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Aims. We constrain the space density and properties of massive galaxy candidates at redshifts of z > 3.5 in the Great Observatories 
Origin Deep Survey North (GOODS-N) field. By selecting sources in the Spitzer + IRAC bands, a highly stellar mass-complete 
sample is assembled, including massive galaxies which are very faint in the optical/near-IR bands that would be missed by samples 
selected at shorter wavelengths. 

Methods. The z > 3.5 sample was selected down to m AB = 23 mag at 4.5 fim using photometric redshifts that have been obtained by 
fitting the galaxies spectral energy distribution at optical, near-IR bands and IRAC bands. We also require that the brightest band (in 
AB scale) in which candidates are detected is the IRAC 8.0 fim band in order to ensure that the near-IR 1.6 fim (rest-frame) peak is 
' falling in or beyond this band. 

Results. We found 53 z > 3.5 candidates, with masses in the range of M* ~ 10 10 - 10" M Q . At least -81% of these galaxies 
are missed by traditional Lyman Break selection methods based on ultraviolet light. Spitzer + MIPS emission is detected for 60% 
of the sample of z > 3.5 galaxy candidates. Although in some cases this might suggest a residual contamination from lower redshift 
star-forming galaxies or Active Galactic Nuclei, 37% of these objects are also detected in the sub-mm/mm bands in recent SCUBA, 
AzTEC and MAMBO surveys, and have properties fully consistent with vigorous starburst galaxies at z > 3.5. The comoving number 
density of galaxies with stellar masses > 5 x 10 10 M Q (a reasonable stellar mass completeness limit for our sample) is 2.6 x 10~ 5 Mpc~ 3 
(using the volume within 3.5 < z < 5), and the corresponding stellar mass density is ~ (2.9 + 1.5) x 10 6 M G Mpc~ 3 , or about 3% of 
the local density above the same stellar mass limit. For the sub-sample of MlPS-undetected galaxies, we find a number density of 
~ 0.97 x 10 5 Mpc 3 and a stellar mass density of ~ (1.15 ± 0.7) x 10 6 M Q Mpc 3 . Even in the unlikely case that these are all truly 
quiescent galaxies, this would imply an increase in the space density of passive galaxies by a factor of ~ 15 from z ~ 4 to z = 2 , and 
Q ■ by ~ 100 to z = 0. 

Key words. Cosmology: observation, - Galaxies: formation - Galaxies: evolution- Infrared: galaxies 
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1. Introduction break galaxies (LBGs) at z — 3, has proven successful also 

for identifying blue star-forming galaxies at higher redshifts 

^ . The question of the galaxy mass assembly at early cosmo- (z ^ 4 _ 6 Steidel et al . 1999 . Dickinson 1998; Steidel et al. 

logical epochs is still open. In a ACDM Universe, the stan- 2003; Giavalisco et al 2004b; McLure et aL 2006; Yan et al. 

dard formation scenario predicts hierarchical growth of struc- 2006) 
tures. Local galaxies would have been formed through repeated 

merging events. The first semi-analytic galaxy-formation models A substantial amount of spectroscopic confirmations is avail- 

(Kauffmann & Chariot 1998) predicted that the assembly of the able for LBGs at z ~ 5 - 6 (Bunker et al. 2003; Stanway et al. 

most massive systems would only have been completed in the 2004; Dickinson et al. 2004; Eyles et al. 2007; Stark et al. 2007; 

most recent epochs (z < 1). Nevertheless, in recent years a sub- Vanzella et al. 2008). This allowed early estimates of the stellar 

stantial population of massive galaxies has been found beyond mass density at those redshifts. For instance Yan et al. (2006) 

z ~ 1. Up to z - 2.5 the most massive objects spectroscopically measured the comoving stellar mass density at z ~ 6 in the total 

confirmed are similar to the local Early Type Galaxies (ETGs): GOODS field and found it in good agreement with the simula- 

old (age of ~1 Gyr), passively evolving, and with M* > 10 11 M tions of Night et al. (2006). The stellar mass density found in 

(Cimatti et al. 2004; McCarthy et al. 2004; Daddi et al. 2005; the GOODS-S field by Eyles et al. (2007) for a sample of i- 

Saracco et al. 2005; Kriek et al. 2006). dropout LBGs in the same redshift range is also consistent with 

At higher redshifts the availability of multi-wavelengths data the Yan et al. (2006) results. However it is larger by a factor of 

from the new generations of deep surveys has allowed searches 4 (~ 2, if the different IMFs are taken into account) with re- 

for massive galaxy candidates up to z ^ 5 - 6.5. The dropout- spect to the Bower et al. (2006) semi-analytical model predic- 

technique pioneered by Steidel et al. (1996) to find Lyman- tions. The general picture is, however, still controversial and fur- 

ther studies are necessary to better constrain galaxy formation 

Send offprint requests to: aaa models. In fact, the Ly man-break technique can only be used 
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to select actively star-forming galaxies having sufficiently lumi- 
nous UV emission to allow their detection and color selection. 
Galaxies, and especially massive ones, could be faint or unde- 
tected in the optical and near-infrared bands. Such sources could 
be either old systems with small quantities of dust and low lev- 
els of star-formation rate, or dust reddened starburst galaxies. If 
such z > 3.5 objects really exist, they would be too faint to be 
studied spectroscopically with the instrumentation that is cur- 
rently available, and the only possible approach that can be used 
to constrain their number density and their properties remains 
a photometric one, based on the Spectral Energy Distribution 
(SED) fitting analysis. 

The most popular criterion used so far to select high-z 
galaxies with red rest-frame colours is based on the redshifted 
4000 A/Balmer break, i.e. the typical features of galaxies with 
evolved stellar populations. It was pioneered by Franx et al. 
(2003) to pick out 'distant red galaxies' (DRGs) at z > 2 
(J s - K s > 2.3, AB system). Recently it was extended by 
Brammer & van Dokkum (2007) to identify DRGs up to z ~ 3.7 
(H - K > 0.9, AB system). By comparing the properties and the 
space density of DRGs at z ~ 3.7 and z ~ 2.4 in the same field, 
they found a stellar mass density of about a factor of 5 lower in 
the higher redshift bin. 

In the recent years the availability of Spitzer data has also 
enabled us to search for galaxies with red rest-frame optical col- 
ors at high redshifts, and to better constrain their stellar masses 
(Fontana et al. 2006; Dunlop et al. 2007; Rodighiero et al. 2007; 
Wiklind et al. 2008; Eyles et al. 2005). Some of the objects iden- 
tified may be among the most massive stellar systems found at 
z ~ 4 - 5. In particular, Rodighiero et al. (2007) extracted a 
Spitzer+IRAC selected sample in the GOODS-South field, lim- 
ited to galaxies undetected in the optical and close to the detec- 
tion limit in the TT-band. Their criterion is complementary with 
respect to those adopted by previous studies (i.e. the K < 23.5 
selection used by Dunlop et al. 2007) and is aimed at identi- 
fying high-z massive galaxies that are missed by conventional 
selection techniques based on optical and near-IR observations. 
They found a potential population of optically obscured massive 
galaxies at z > 4. In the same field Wiklind et al. (2008) found 
11 ^-selected z > 5 massive and evolved galaxy candidates. 
They selected these objects by applying color criteria aimed at 
the identification of the 4000 A/Balmer break between the K and 
3.6 pm passbands, in combination with SED fitting. 

It should be noted that in both of these works roughly half of 
the high-z massive candidates are detected at 24 pm. The 24 pm 
emission might be due to obscured AGN activity at high redshift, 
but also from Polycyclic Aromatic Hydrocarbon (PAH) emis- 
sion in dusty star forming galaxies at lower redshifts. In fact, 
one has to take into account the degeneracy between reddening 
and redshift while selecting high redshift sources based on pho- 
tometry. Recent studies of sub-mm galaxies (SMGs) have shown 
that moderate 24 /im emission can arise from massive starburst 
galaxies at high redshift (z > 3.5, Daddi et al. 2008; Greve 
et al. 2008; Perera et al. 2008). Hence one of the main difficul- 
ties in selecting bona fide high-z galaxies through photometric 
redshifts is to efficiently remove the contamination from dust- 
reddened galaxies at lower redshift. As an example, one of the 
most debated cases is the galaxy HUDF-JD2 in the Hubble Ultra 
Deep Field. Mobasher et al. (2005) suggested that such a galaxy 
is a massive (M 6 x 10 11 M ), low-reddened, old galaxy (with 
an age of ~ 10 8 yrs and a Zf orm > 9) at z - 6.5. However, other 
studies (e.g. Dunlop et al. 2007; Rodighiero et al. 2007; Chary 
et al. 2007) suggest that HUDF-JD2 is instead a lower redshift 



galaxy with very strong dust reddening (z ^ 2.2 and A v - 3.8, 
Dunlop et al. 2007). 

The principal aim of this paper is to study the properties 
and space density of massive galaxies at z > 3.5 in the Great 
Observatory Origin Deep Survey-North (GOODS-N). This sur- 
vey is especially suitable for our purpose because of the quality 
and depth of the available multi-wavelength data. In particular it 
has the advantage of a better coverage at sub-millimetric wave- 
lengths with respect to GOODS-S. In addition, while similar 
searches have been carried out in the GOODS-South field, none 
has been performed so far in the northern hemisphere field of 
the GOODS Survey. We chose to select the sample in the IRAC 
channel 2 (4.5 pm) for two reasons. First, to pick out the most 
massive systems. At z > 3.5 the 4.5 pm band samples the rest- 
frame near-IR emission, which is well correlated to the galaxy 
stellar mass (Kauffmann & Chariot 1998). Second, the 4.5 pm 
band selection allows us to have a more complete sample with 
respect to an optical or near-IR selection. This is one of the main 
differences with the previous works: we attempted to construct 
a sample as complete as possible, including objects undetected 
in optical and near-IR bands, and hence missed by the optical 
or K selection criteria used so far in the literature (Drory et al. 
2004; Fontana et al. 2006; Dunlop et al. 2007; Wiklind et al. 
2008; Brammer & van Dokkum 2007). As already mentioned, 
Rodighiero et al. (2007) also used an IRAC band (3.6 pm) se- 
lection, but their additional conditions about the lack of optical 
emission and the faint near-IR detection (K > 23, AB system) 
make their sample incomplete by construction. 

The structure of the paper is as follows. In Sect. 2 we present 
the optical, near-IR and IR data available for the GOODS-N 
field. In Sect. 3 we describe the IRAC selection and the photom- 
etry of the sample. Sect. 4 describes the results from the SED 
fitting analysis with the measurement of photometric redshift 
and the stellar mass estimates. In Sect. 5 we discuss the selec- 
tion of z > 3.5 galaxy candidates. In the same section we dis- 
cuss the detections of the sample in the MIPS-24 pm and in the 
sub-mm/mm. In Sect. 6 we study the effects of incompleteness 
and the bias against low mass galaxies in our magnitude-limited 
sample. We give estimates for the comoving number density and 
stellar mass density for the full sample, and for the sub-sample 
of MIPS -undetected galaxies. Sect. 7 presents the summary and 
conclusions. 

Throughout this paper, we adopt the following cosmological 
parameters Ho = 70km s _1 Mpc -1 , Qa = 0.7 and Cl m = 0.3, and 
all magnitudes are given in the AB system. 

2. The GOODS-N data. 

This work is based on the analysis of the Great Observatory 
Origins Deep Survey - North (GOODS-N) data. The GOODS- 
N is a multi-wavelength deep survey centered in the Hubble 
Deep Field North (HDFN 12 /, 36 m 55 ? ,+62°14'15" Giavalisco 
et al. 2004b) with size of ~ 10' x 16'. The field is part of 
the GOODS Spitzer Legacy Program (PI: M. Dickinson) and 
the Hubble Space Telescope (HST) Treasury Program (PI: M. 
Giavalisco), divided in two fields (GOODS-N and GOODS-S), 
from the northern and southern hemispheres (see Dickinson et al. 
2003). The GOODS-N field observations span from the X-rays 
to the sub-mm wavelengths. In particular, for the analysis pre- 
sented in this paper, we make use of the following data sets: 

- U band imaging, taken at the Kitt Peak National Observatory 
(KPNO) 4m telescope, from the Hawaii Hubble Deep Field 
North (H-HDFN) survey datasets (Capak et al. 2004). 
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- Optical imaging from the Hubble Space Telescope (HST) 
with the Advanced Camera for Survey (ACS) in the 
F435W (B), F606W (V), F775W (i), F850LP (z) pass-bands 
(Giavalisco et al. 2004b). 

- Near-Infrared imaging, in the J, H and Ks bands, from 
the KPNO-4m telescope (M. Dickinson private communica- 
tion). 

- Infrared observations from the Spitzer Space Telescope 
(Werner et al. 2004) + Infra-Red Array Camera (IRAC; 
Fazio et al. 2004) 

- 24 pm observations from the Spitzer Space Telescope and 
Multi-band Imaging Photometer for Spitzer (MIPS; Rieke 
et al. 2004). 

Specifics of the data are given below. 

2.1. KPNO 4m U-band imaging 

The H-HDFN Survey is a deep multi-colour imaging survey of 
0.2 degrees 2 , centered in the Hubble Deep Field North. The data 
consist of deep images in the U, B, V, R, I, and z' bands. Out of 
these we used only the U-band images, due to the availability 
of HST/ ACS imaging from the B through z bands. The U-band 
observations were taken in March 2002, using the KPNO-4m 
telescope with the MOSAIC prime focus camera (Muller et al. 
1998), covering a field of view of 36' x36' with a seeing of 1'.'26. 
We analyzed the part overlapping with the GOODS-N field (an 
area of about 165 arcmin 2 around the HDF-N). The U image 
reaches a 5<x sensitivity limit of 27. 1 mag. 

2.2. ACS+HST optical imaging 

The GOODS-N optical data from the Advanced Camera for 
Surveys (ACS), on-board HST, are in the F435W (B), F606W 
(V), F775W (i), F850LP (z) pass-bands (vl.O data products). 
The 5cr point source sensitivity limits are 28.5, 28.8, 28.1, 27.6 
respectively, scaled from values reported in Giavalisco et al. 
(2004b) to the full exposure time of ACS vl.O images. The 
HST/ ACS high-resolution data is organized in 17 images, called 
"sections". Each section is an image 8 192 x 8 192 pixels in size, 
drizzled from the native ACS pixel scale (0705/pixel) to an im- 
age subsampled at 0703/pixel. The huge size (~ 40 000 x 40 000 
pixels) of the combined HST images in the field made source ex- 
traction an extremely intensive computing task. Therefore we re- 
binned the final image using an 8 x 8 pixel kernel, requiring flux 
conservation. The new re-binned mosaic obtained in this way 
has a size of 5 120 x 5 120 pixels, and a pixel scale of 0724/pixel. 

2.3. KPNO 4m FLAMINGOS near-infrared imaging 

We complemented the public data with near-IR imaging taken at 
the KPNO-4m telescope with the Florida Multi-object Imaging 
Near-IR Grism Observational Spectrometer (FLAMINGOS; ob- 
servations by Mark Dickinson et al.). The FLAMINGOS images 
in the J, H and Ks pass-bands, are mosaics of data taken at dif- 
ferent pointings and orientations to cover the whole GOODS- 
N. The average exposure times per band in the main part of 
the GOODS-N area are 20250 s, 15770 s and 35440 s at J-, 
H- and Ks, respectively. The three images reach a 5 cr point 
source sensitivities of 24.03, 23.77, 23.81 mag. The pixel-scale 
is 07316/pixel, and the seeing 1727, 1'.'2 and 172, for the J, H, 
and Ks pass-bands, respectively. 

We performed some photometric checks before measur- 
ing source fluxes. We verified the zero-points reliability of 



FLAMINGOS images, comparing the bright and unsaturated 
star magnitudes with stars from the 2MASS catalogue. 

2.4. Spitzer imaging at 3.6 to 24 pm 

The infrared imaging is from Spitzer Space Telescope + IRAC 
between 3.6 and 8.0 pm, and from Spitzer + MIPS at 24 pm. 

The four IRAC channels, centered on 3.6 pm, 4.5 pm, 
5.8 pm, 8.0 pm, have formal 5 cr limits for isolated point sources 
of 26.1, 25.5, 23.5, 23.4 respectively. The mean FWHM of a 
point source for the IRAC-bands are V!66, 1'.'72, 1788 and 1798, 
for channels 1, 2, 3 and 4 respectively. The initial image pixel- 
scale of 1722 is reduced to 076/pixel, after the dithering and driz- 
zling process. 

The MIPS public dataset includes calibrated maps and a cat- 
alogue of 24 pm sources with flux densities S 24 > 80 piy. 
The PSF was generated from isolated sources in the image, and 
renormalized based on the aperture correction published in the 
MIPS data Handbook (v2.1, section 3.7.5, table 3.12). Then an 
independent PSF algorithm was run to extend the 24 pm sample 
to fainter sources (see Rodighiero et al. 2006). By this procedure 
the detection threshold was extended to S24 > 35 piy. Since the 
24 pm information can help to better constrain the SFR (Chary & 
Elbaz 2001) and the nature of high-redshift sources, we also ex- 
ploited a deeper 24 pm catalogue, reaching a 5cr limit of 20 yuJy 
(Chary et al. 2008 in preparation, Daddi et al. 2007). 

3. Sample selection 

In order to have an accurate measurement of the object fluxes in 
all bands, we performed aperture photometry in each band us- 
ing the SExtractor software (Bertin & Arnouts 1996). The initial 
catalogue was built by selecting the sources in the 4.5 pm pub- 
lic image obtained with Spitzer Space Telescope + IRAC, down 
to a limiting magnitude of 1714,5 = 23.0 (2.3 ^Jy). Simulations 
indicate that a sample selected in this way is complete at the 
80% level at 1714s ~ 23.0 (M. Dickinson private communication), 
where 20% of the sources are lost due to blending. 

The infrared selection allows us to search for the presence of 
massive galaxies up to very high redshifts. In fact for galaxies at 
3.5< z <7 the rest-frame light emitted at A > 0.6 pm is redshifted 
into the IRAC 4.5/im band. Therefore, the 4.5pm selection is 
directly sensitive to stellar mass, rather than to the ongoing star 
formation activity. From a technical point of view, the choice of 
selecting the sample at 4.5 pm represents the best compromise 
amongst the IRAC bands in terms of image quality, blending 
problems, and sensitivity. 

In order to detect galaxies in the 4.5/im IRAC band and per- 
form the photometry we used SExtractor with the following set 
of parameters. The detection limit was set at ~ 1 cr above the sky 
background. A Gaussian filter was used to improve the detection 
of faint extended objects. Due to the high crowding of IRAC im- 
ages, and the large PSF, the minimum contrast parameter was set 
to a small value of 5 x 10~ 9 , to help with source deblending. 

We measured fluxes in 470 diameter apertures, as this size 
allows us to minimize uncertainties in the photometry in the 
Spitzer +IRAC bands (as also suggested by the SWIRE team 1 ). 
In the 4.5 pm IRAC filter we applied an aperture correction of 
0.25 mag that we independently computed from point-source ob- 
jects, by measuring the total flux in 870 diameter apertures. 

The final IRAC selected catalogue contains 4 142 objects 
(m 4 . 5 < 23.0). 

1 http://data.spitzer.caltech.edu/popular/swire/20050603_enhanced_vl/ 
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3. 1 . Multi-band Photometry 

We searched for the counterparts of the 4 142 IRAC-selected ob- 
jects in the other datasets. We detected the sources and mea- 
sured the photometry independently in each dataset and associ- 
ated counterparts to the 4.5 fim IRAC detected galaxies using a 
search radius of 1" 0. 

For the U-band we measured fluxes 4'.'0 (diameter) apertures. 
We applied an aperture-correction of 0. 1 mag, computed by con- 
sidering the total flux in 8'.'0 diameter aperture. For undetected 
sources we used the 3 cr flux upper limit of 26.7 mag. This value 
was computed by measuring the standard deviation of the sky 
signal at random position over the field. 

For the HST/ACS bands circular apertures with diameters 
of 2'.'0 were used. The following aperture-corrections for point- 
sources were found by measuring the total flux in 5'.'0 diameter- 
aperture, for the B, V, i, z bands: 0.050, 0.040, 0.030, 0.045 mag, 
respectively. For undetected objects we obtained 3 cr upper lim- 
its in the same way as described above for the U-band. We mea- 
sured upper limits to be 26.3, 26.9, 26.3, 25.5, for the B, V, i and 
z bands respectively. 

The near-IR magnitudes were measured in 4'.'0 diameter 
apertures, and corrected to total (within 870 diameter) by 0.1, 
0. 14, 0. 15 mag for the J-, H- and Ks- band, respectively. The 3 cr 
detection limits are 23.3, 23.0 and 23.1, for the J-, H- and Ks- 
band. 

In the other IRAC channels we used the same aperture di- 
ameters as in IRAC-4.5 fim (4'.'0), and the same approach to 
compute aperture corrections. The final IRAC magnitudes were 
obtained by applying aperture corrections of 0.23, 0.25, 0.40, 
0.44 mag for 3.6, 4.5, 5.8 and 8.0 pm bands, respectively. 

To check the reliability of our photometry we compared 
the colors of stars with those expected based on the models of 
Lejeune et al. (1997) in color-color diagrams involving one band 
from each one of the three datasets (ACS/HST optical, Near- 
Infrared/FLAMINGOS and infrared/IRAC-5/?/fzer). The case of 
(b - J) vs {J - m-3.6) is shown in Fig. 1. The star sequence is 
clearly visible on the left side of the panel. As it appear from the 
figure we found a very good agreement between the observed 
colours (green circles) and the expected ones (red asterisks from 
the "corrected" templates of Lejeune et al. 1997) for spectro- 
scopically identified stars (public database of the TKRS survey 
Wirth et al. 2004). Therefore, thanks to the very good concor- 
dance between optical, near-IR and IR datasets we did not have 
to apply further photometric corrections to our data magnitudes. 

The photometric errors given by SExtractor software are 
generally underestimated, because the software does not take 
into account the possible correlation between neighboring pixels 
that can result from image processing. Hence we applied cor- 
rections to the photometric errors in each data-set. For the four 
IRAC channels we added 0. 1 mag quadratic to SExtractor errors. 
This is a conservative choice to account for uncertainties on the 
photometric zeropoints and on the aperture corrections, as also 
suggested in Maraston et al. (2006). In the other filters we mea- 
sured the signal in a number of random sky positions in the field, 
within the same aperture used to perform galaxy photometry in 
each band. Comparing the standard deviation of these measure- 
ments to the actual errors in the photometry of faint galaxies 
inferred by SExtractor we derived scaling factors to apply to the 
formal errors reported in the SExtractor catalogues. 
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Fig. 1. The color-color diagram (b - J) vs (J - m^s). The star sequence 
is clearly visible on the left side of the diagram. We found a very good 
agreement between the observed colours (green circles) and the ex- 
pected ones (red asterisks from the "empirically corrected" templates 
of Lejeune et al. 1997) for spectroscopically identified stars (public 
database of the TKRS survey Wirth et al. 2004). 

4. Spectral Energy Distribution (SED) fitting 
analysis 

4.1. Photometric redshifts estimate 

To estimate the photometric redshifts of our sample of galax- 
ies we used the Hyperz code (Bolzonella et al. 2000). This code 
compares the multi-wavelength photometry of each source to a 
database of theoretical and/or empirical templates at different 
redshifts. The best fit solution, for a given set of templates, is 
derived through a standard^ 2 minimization. 

For the determination of the photometric redshifts we ex- 
cluded the 24 fj.m band, because the templates of stellar popu- 
lations that we used do not include dust emission. For a similar 
reason we decided to omit also the 8 pm band, which includes 
thermal dust and PAH (i.e. non-stellar) emission for the lower 
redshift galaxies and can be contaminated by star formation or 
AGN also at higher redshift (Daddi et al. 2007). We used the 
template library of the observed SED constructed by Coleman 
et al. (1980, hereafter CWW), extended in the UV and IR re- 
gions by means of Bruzual & Chariot (2003) synthetic spectra, 
as explained in the Hyperz User's manual 2 . This set of templates 
includes four types of spectra with characteristics correspond- 
ing to the different local types of galaxies: Elliptical/Lenticular 
(E/S0), Spiral (Sbc and Scd) and Irregular (Im). The CWW li- 
brary was complemented with a spectrum of a Young Star-Burst 
(SB) galaxy, with constant star formation rate (SFR) and age of 
0. 1 Gyr, from the Bruzual & Chariot (2003) models. 

We took into account the dust extinction by applying the 
Calzetti et al. (2000) law. The extinction parameter (Ay) was al- 
lowed to vary from to 0.8 in steps of d(Ay) = 0.1. Only a small 
range of Ay was allowed since the CWW templates include al- 

2 http://webast.ast.obs-mip.fr/hyperz/ 
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ready intrinsically reddened SEDs. The permitted redshift range 
span the large interval of values of < z < 10, with a step of 
d(z) = 0.07. 

We used the prescription of Madau (1995) to represent the 
average Lyman-Q' forest opacity as a function of wavelength and 
redshift. 

4.2. Photometric redshift reliability 

We checked the reliability of our photometric redshifts for 836 
objects of the GOODS-N field having high quality spectroscopic 
redshifts. Unobscured AGN were excluded from the list by 
means of the X ray information from the publicly available X-ray 
catalogue of the Chandra Deep Field-North Survey Alexander 
et al. (2003). The largest number of spectroscopic redshifts 
(721) are from the public database of the Team Keck Treasury 
Redshift Survey (TKRS; Wirth et al. 2004), and are accurate 
to 100 kms~ l . The median redshift of this survey is rather low 
(z = 0.65). In order to increase the number of available objects 
with the highest spectroscopic redshifts, we complemented the 
TKRS database with the galaxies found by Reddy et al. (2006) 
in the redshift interval of 2.0 < z < 3.5 (108). We also used 
some objects (5) up to z — 5 from the NASA/IPAC Extragalactic 
Database 3 (NED). 
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Fig. 2. Comparison between photometric z p i,ot and spectroscopic z spe c 
redshifts for 836 objects of our sample. The z spec -Zphot relation has a 
combined mean offset of (z spe c - Zpfc()/(1 + Zspec) — 0.004 with an rms 
scatter of cr [(z, spec - z p i,ot)/(l + z spec )\ = 0.09. For the flagged objects 
(red open diamonds), appearing as catastrophic outliers we suggest the 
lower redshift solution, as detailed in the text. 

The comparison between photometric (z p hor) and spectro- 
scopic (Zspec) redshifts is shown in Fig. 2. After clipping 5 cr out- 
liers we found a combined mean offset of (z sp ec - ZphotMi^ + 
Zspec) = 0.004 and arms scatter of tr [(z spe c-Zphot)/(I +z spe c)] = 
0.09. This result was obtained by means of an iterative optimiza- 
tion procedure applied to the Hyperz results. It consists in re- 
peating the SED fitting analysis for the same objects, fixing the 
redshifts to the known spectroscopic values. For each filter the 
median offset between empirical and "theoretical" magnitudes 
(the latter deriving from the best fit SEDs) was obtained. The 
largest photometric offsets that we found were for the U -, z- and 
/T-bands (+0.2, +0.09, +0.14, respectively) and | Amag|~ 0.05 

3 http ://nedw w w. ipac.caltech.edu/ 



Table 1. Outliers of our z p i m t - z spec comparison (Fig. 2). Their best 
fit SEDs are shown in Fig. 3. The coordinates and the spectroscopic 
redshifts are from NED database. The last column is the distance to the 
NED position (arcsec). 



id 


RA 


Dec 


Zphot 


Zspec 


dist(") 


#6955 


12:36:40.510 


+62:13:34.90 


1.48 


3.826 


0.27 


#8734 


12:36:37.630 


+62:14:53.70 


0.68 


4.886 


0.16 


#8903 


12:37:21.048 


+62:15:01.73 


0.98 


4.762 


0.41 



for the others. Adding these corrections to the observed magni- 
tudes, the rms of the A z/(l +z sp ec) quantity was slightly reduced. 
The same offsets were finally applied to the whole sample for es- 
timating the photometric redshifts. 

Since the current work is focused on galaxies at z > 
3.5 we have carefully examined three outliers (#6955, #8734 
and #8903) in the photometric-spectroscopic redshift compari- 
son with reported z spec > 3.5. These are noted with red diamonds 
in Fig. 2, and information about them is presented in Tab. 1. In 
Fig. 3 we show the photometric SEDs for these three objects and 
in Fig. 4 we show cutout images from U -band through 24 //m. 
For the first two objects (#6955 and #8734) Dawson et al. (2001) 
report Keck spectra showing isolated emission lines with no con- 
tinuum detection (their quality class 4; see Tab. 2 and 3 of that 
paper), which they identify as weak Ly-a emission. However ob- 
ject #6955 is detected at a significance > 3 <x limit in the U -band, 
ruling out the proposed z = 3.826. More recent Keck/DEIMOS 
spectra from 2005 and 2006 (M. Dickinson & D. Stern, private 
communication) detected an emission line at 9273A which is 
likely to be [OII]3727A at z = 1.488, consistent with our pho- 
tometric redshift estimate for this object. Object #8734 shows a 
red point source in the HST/ACS images, along with faint, bluer, 
extended emission to the southeast. This is very likely to be a 
superposition of a cool galactic star with a faint galaxy. Cowie 
et al. (2004) report that this is a star based on DEIMOS spec- 
troscopy. The faint, extended component is certainly a galaxy, 
although it is clearly detected in the ACS Z?-band image, and 
is therefore unlikely to be at z = 4.887 as reported by Dawson 
et al. (2001). Overall, the photometry at all wavelengths redward 
of the V-band is dominated by the point source component, and 
the photometric SED for this object Fig. 3 clearly resembles that 
of a cool, late-type star (probably an M-dwarf). Finally, for ob- 
ject #8903, Dawson et al. (2001) report z = 4.762. However, 
ACS V - i and i - z colors for this galaxy are not consistent 
with it being a V- dropout Lyman break galaxy, and fall within 
the locus of ordinary foreground galaxies. The ACS color im- 
ages show a faint, red, early-type galaxy, and our photometric 
redshift estimate is z p hot = 0.89, very close to that proposed by 
Capak et al. (2004). Although Cowie et al. (2004) do not report 
details about the spectral features that they observed, it is possi- 
ble that faint [Oil] emission at z = 0.88 was mistaken for faint 
Lya at z — 4.76. 

4.3. Galaxy mass estimate 

We used the code Hyperzmass to estimate the mass for each of 
our sources. The Hyperzmass code works in a similar way to 
Hyperz(see Pozzetti et al. 2007). In this task the redshifts were 
fixed at the derived photometric values, or spectroscopic ones, 
when available. The observed data points were fitted using both 
the synthetic spectral models of Bruzual & Chariot (2003, here- 
after BC03) and Maraston (2005, hereafter MA05). Both tem- 
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Fig. 3. SEDs of the outliers in the z spe c-Zphot space (see objects flagged 
in Fig. 2). Our best fits (black curves) produce photometric redshifts of 
~1.48, 0.68 and 0.98, respectively. We compare these results with the 
best fits (red curves) obtained by fixing the redshifts to the spectroscopic 
values taken from the literature (z = 3.826, z = 4.886, and z = 4.762; 
NED database). Obviously, the black curves fit the observed data better 
than the red ones. See also Tab. 1 



plate sets are composed by models of stellar populations with 
Star Formation Rates (SFRs) exponentially decreasing with time 
(if/(t) = T~ l exp~'l T ). In our analysis we considered models with 
different values of the time-scale r (in unit of Gyr) : 0.1, 0.3, 1, 



U #6955 #8734 #8903 




Fig. 4. Multi-wavelength identification of the three outlier objects in the 
z S p e c-Zphot diagram. The cut-outs have a size of 10" x 10". North is up. 
East is at left. 
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2, 3, 5, 10, 15, 30 and oo. The last model corresponds to a stellar 
population with a constant SFR. 

For the BC03 models we used the Chabrier (2003) Initial 
Mass Function (IMF). Instead, the MA05 models adopt the 
Kroupa (2001) IMF, which is relatively similar to the Chabrier 
IMF. To compare the different results we computed a systematic 
correction to convert the masses obtained with Kroupa IMF to 
Chabrier IMF (log M cha brier = log M Kroup? - 0.04). 

In the SED fitting analysis the extinction parameter, Ay, was 
allowed to span the range of values < Ay < 6, with a step of 
0.1. This large reddening range is required given that synthetic 
stellar population models are not intrinsically reddened as the 
observed ones. We used Ay = 6 as an upper limit for the red- 
dening, since there is no evidence of objects with A y higher than 
~ 5. The galaxy age was left as a free parameter, spanning the 
range from 5x 6 yr to 11.7 x 10 9 yr (and always requiring ages 
less than the age of the Universe at each redshift). To minimize 
the number of free parameters we limited the fits to solar metal- 
licity. 

5. Massive galaxies at z > 3.5 

As previously discussed, evidence for populations of massive 
galaxies at z > 3 has been found by several studies. However, 
results are not yet conclusive since the number of objects is low, 
and frequently biased against objects that are faint or undetected 
in optical and near-IR bands. Moreover for most of these objects 
no spectroscopic redshift information is available. Also, the pos- 
sible contamination by lower redshift interlopers is a source of 
uncertainty. To shed some light on the nature and the number 
density of the high redshift population of massive galaxies we 
selected a final sample of z p hot > 3.5 candidates, by applying 
the following selection criteria (summarized in Table 2). Among 
the initial 4 142 IRAC-4.5 pm selected galaxies we found 190 
galaxies with photometric redshifts of z p hot > 3.5. 

All such objects were carefully inspected for possible blend- 
ing issues. One of the main problems in using Spitzer data in 
crowded fields is, in fact, the blending between nearby sources. 
Neighboring objects that can be well separated in optical and 
near-IR images, might happen to be blended in the IRAC bands, 
due to the limited spatial resolution of the Spitzer telescope. 
Hence, in order to build a bona-fide sample of high-z galaxies 
we had to take into account that a fraction of the IRAC-selected 
objects could be affected by blending. We decided to visually in- 
spect all the z > 3.5 objects, and we excluded from the sample 
any sources that the SExtractor software was not able to deblend. 
To identify them we visually inspected examined the SExtractor 
"aperture image" in each filter. This is a 'check-image' given as 
output by the software, where the circular apertures are shown 
for all the objects detected above the threshold level. We classi- 
fied an object as blended if the IRAC position did not correspond 
to a single source in the higher resolution optical and/or near- 
IR images (e.g. objects that are well separated in the optical or 
near-IR images may be blended in the IRAC images due to their 
poorer angular resolution). About the 25% of the pre-selected 
Zphot ^ 3.5 objects were flagged as blended and therefore dis- 
carded. 

By matching our catalogue with the publicly available X-ray 
catalogue of the Chandra Deep Field-North Survey 3 , we found 
14 objects that were detected in the hard X-ray band, and one 
with soft X-ray flux. We excluded these 15 sources from the final 
sample since X-ray emission indicates the presence of AGNs in 
those galaxies, that can alter the observed flux and induce an 
overestimate of the stellar masses. 



Finally we introduced a further criterion, based on IRAC 
colours. We considered as genuine z > 3.5 candidates only those 
sources showing an SEDs peak at the 8.0 pm observed frame 
(1-2 pm rest-frame for z > 3.5 galaxies) in the AB scale. The 
reason is that the near-infrared emission of all but extremely 
young galaxies show a distinct peak at 1 .6 pm in the rest-frame, 
that is redshifted in the 8.0 pm IRAC band for z > 3.5 sources. 
This peak corresponds to the combination of the Planck spec- 
tral peak of cool stars and the effects of a minimum in the H 
opacity in the stellar atmospheres. It has been used as photomet- 
ric redshift indicator in several previous studies (e.g. Sawicki 
2002; Berta et al. 2007). This criterion was applied because at 
this step we expect that our high-z sample will be contaminated 
by a substantial fraction of lower-redshift galaxies, because of 
the intrinsic uncertainty of the photometric redshift estimate. 
The degeneracy between reddening and redshift could introduce 
lower-redshift contaminants in the final high-z sample (i.e. dust 
reddened galaxies at lower redshift can be mistaken for higher-z 
old stellar populations). Given that galaxies at z > 3.5 represent 
only a minority (~ 4%) of the galaxy sample selected at 24/im, 
while most of the objects are at z < 2, from a statistical point 
of view we could expect a more substantial contamination by 
lower-z sources entering our z > 3.5 redshift selection, largely 
overbalancing the number of genuine z > 3.5 galaxies lost as 
due to errors in their photometric redshifts. In addition, as al- 
ready pointed out, the 8 pm data-point was excluded from the 
SED fitting analysis because for objects at z < 1.8 it could sam- 
ple the PAH emission lines for dusty galaxies, and the stellar 
population models that we used do not include dust emission. In 
this way, by excluding the 8 pm band from the fitting procedure, 
we may fail to distinguish galaxies at z ~ 2 - 3 where the 1 .6 pm 
peak lies in the 5.8 pm band from those at z > 3.5. Requiring the 
SED to rise from 5.8 to 8 pm can help to exclude these possible 
lower-redshift contaminants. 

The parameterization of the colour criterion we used is the 
following: 

(m 8 .o < m 5 . 8 + 0.1) n (m 8 .o < m 4 . 5 + 0.1)n (1) 
n(m 8 . < m 3 . 6 + 0.1) 

the 0. 1 mag terms are introduced to account for the median un- 
certainty in the IRAC photometry of ~ 0. 1 . 

In the two panels of Fig. 5 we show the expected ms.s - ffig.o 
colour as a function of redshift for different stellar population 
models. The top panel shows Simple Stellar Populations (SSPs) 
with three different ages from 0.1 to 1 Gyr, and the bottom 
panel shows a constant star formation rate model with fixed 
age of 0.1 Gyr and reddening from E B -v =0 to 0.8. The tem- 
plates shown are from the Maraston stellar population models 
(Maraston 2005). These models, however, are not supposed to 
be used for A > 2.5 pm where they display a strong, unphys- 
ical, discontinuity. This discontinuity affects the IRAC colours 
at < z < 2.7. Hence we tested the reliability of this colour 
criterion also with BC03 templates, that are not affected by this 
problem. We verified that the expected IR colours for BC03 stel- 
lar population models at z < 2.7 are equally bluer than the colour 
cuts defined by the equation 1 . 

From the upper panel of Fig. 5 it is clear that the 1.6 pm 
peak criterion may lead to exclude some blue and unreddened 
(E B -v = 0) galaxy at z > 3.5 with age <0.1 Gyr. This is be- 
cause galaxies dominated by young stellar populations do not 
show the 1.6 pm peak in their SEDs (see Sawicki 2002). In fact, 
due to this criterion we excluded from the final sample a spectro- 
scopic confirmed LBG (object #5455, from Spinrad et al. 1998, 
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with ispec = 5.34 and z p hot = 5.18), for which the best fit re- 
sult predicted age ~ 7 Myr, and mass ~ 10 10 M . Nonetheless 
we decided to apply this criterion as a conservative way to ob- 
tain a reliable sample of the most massive galaxies at z > 3.5. 
The objects missed in this way are very few (see next section), 
and are also the least massive ones. Their inclusion would not 
strongly affect our estimate of the comoving stellar mass density 
(for more details see Sect. 5.2, 6.1 and 6.2). 



Table 2. Sample selection criteria. In the left and right columns the sub- 
sequent steps to select the final sample and the corresponding number 
of objects are shown, respectively. 



Selection Criteria 


Number of Objects 


m 4 .5 < 23 


4142 


Zphot > 3.5 


140 


No hard X-ray emission 


125 


IRAC/channel 4 peak 


53 



: SSP, MA05 moc 


els : 
— — " 1 Gyrs - 


^ 1- — ' — 0.3 Gyrs: 




OlfGyrs! 



o - 



I COST SFR, MA05 models I 


I age = 0. 1 Gyr 


E_b-v = 0.8 : 




- _--E---b--^ = 0,.5_: 




€^b = v = 0.3_- 




' ' E_b_-v = - 


: \ v/ 





Fig. 5. The two panels show the m 5 8 - m 8 .o colour in function of redshift 
for the Maraston stellar population models (Maraston et al. 2006). Top 
panel: Single-burst Simple Stellar Populations (SSPs) with three differ- 
ent ages (0. 1 , 0.5 and 1 Gyr). Bottom panel: 0. 1 Gyr stellar populations 
with constant star formation rate (SFR) and four different extinction 
values (E B _ V =0.0, 0.3, 0.5, 0.8) 



This criterion is very efficient in excluding SSPs galaxies 
with redshifts lower than ~ 2.5. On the other hand, we should 
point out that it cannot completely remove the contamination by 
interlopers at 2.5 < z < 3.5, as SSPs models with ages > 1 Gyr 
and 3 < z < 3.5, and dusty star forming galaxies at 2.5 < z < 3.5 
can still satisfy this criterion (see the bottom panel of the same 
figure). 

Finally it should be noted that the adopted colour criterion 
would be met also by the so-called TRAC power-law sources'. 
These objects are generally selected by requiring a rising SED 
over the four IRAC bands, and they are believed to be generally 
z < 2.5 AGNs with a hot dust component swamping the stellar 
emission (Alonso-Herrero et al. 2006; Donley et al. 2007). Since 
Donley et al. (2007) selected a sample of (62) IRAC power- 
law galaxies in the Chandra Deep Field Nord (CDFN), whose 
27 are in the GOODS-N field, we investigate the overlap with 
their sample. We found 10 sources in our 4.5 /im IRAC-selected 
(m 4 .5 < 23) sample (of 4142 objects) in common with the TRAC 
power-law' sample. All these sources but one have the infrared 
part of the SED peaked at 8.0 fim, but most of them had been ex- 
cluded from our z > 3.5 final sample because of their lower pho- 
tometric redshifts. Only one object of the Donley et al. (2007) 
sample fulfill all of our selection criteria. This is object #14793 
in our sample (see Tab. 4 and Fig. 14), corresponding to the ob- 
ject CDFN: [DRP2007] 19080 in Donley et al. This galaxy has 
no X-ray emission (< 4.41 x 10~ 16 erg s _1 crrr 2 ) and a rather 
high 24 /im flux (111 ±3.5juJy). The high redshift solution is still 
consistent with these properties. 

As result we obtained a final sample of 53 galaxies (~ 43% 
of the pre-selected z > 3.5 sample), that we consider to be valid 
candidates at z > 3.5. About the 12% of the 125 pre-selected 
high-z candidates have been discarded due to the low S/N ratio in 
the 8.0 //m-band, but we cannot exclude that some of them could 
be at high redshift. The other 45% is likely mainly composed 
either by younger systems (age < 0. 1 Gyr) at the same redshift, 
or by lower redshift z < 3.5 contaminants. We will show in the 
following that the former are only a minority (see Sect. 5.2). In 
Table 6, the multi-band photometry of the final sample is shown. 

As a further test on the quality of this sample we checked 
how many candidates satisfy the empirical blending criterion 
tested on extensive simulations by the GOODS team and used 
by Daddi et al. (2007) and Dickinson et al. (2008 in prep.). This 
is based on the measurements of the angular separation (AO) be- 
tween IRAC positions and the /T-band counterparts (we used 
also optical coordinates for ^-undetected sources). According 
this criterion, if a galaxy has AO > ff!5, it is likely contaminated 
by blending in the IRAC bands. However we should point out 
that this criterion proved to be successful on objects detected at 
the 5-10 cr in both of the K and IRAC bands. For sources with 
fainter optical and near-IR detections the shift of the centroid 
position throughout the different bands could be larger due to 
higher noise. Hence a larger AO could derive from coordinate 
fluctuations rather than from blending contamination. Moreover 
due to the extremely red colours of our candidate high-z galax- 
ies, this criterion is not extensively applicable to our full sam- 
ple being half of it both optical and near-IR undetected. Out 
of the sources with at least a reliable optical or near-IR de- 
tection (>3-5 cr) we verified that the majority (18/25) did pass 
this criterion. On the other hand we found 7 galaxies (#1098, 
#2172, #2796, #13857,#15268,#15541,#15761)with A<9 > 075. 
However for these objects the lack of visual evidence of blend- 
ing led us to conclude that the quite large A6 separation, is a S/N 
issue. Moreover for two of these galaxies the confirmed spectro- 
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scopic redshifts (#13857, #15541, see Sect. 5.4) prove that the 
measured IRAC photometry leads to reliable SED fitting results. 

5.1. Comparison between the BC03, MA05 models 

The SED fitting analysis with BC03 and MA05 templates gives 
somewhat discordant results for both galaxy masses and ages 
estimate. This is an issue recently debated by several authors 
(e.g. Maraston et al. 2006; Kannappan & Gawiser 2007; Bruzual 
A 2007), the major source of this discrepancy being the different 
treatment of the Thermally Pulsing Asymptotic Giant Branch 
(TP-AGB) phase of stellar evolution in the population models. 

In the MA05 models stars in the TP-AGB phase contribute a 
dominant source of bolometric and near-infrared light for stellar 
populations in the age range 0.2 to 2 Gyr, much larger than is 
found in the BC03 models. For those models the TP-AGB phase 
has been calibrated with local stellar populations. Another differ- 
ence is in the treatment of convective overshooting. In their work 
Maraston et al. (2006) tested both the BC03 and MA05 models 
on a sample of high-z galaxies, with reliable spectroscopic red- 
shifts (1.4 < z < 2.5) in the Hubble Ultra Deep Field (HUDF). 
They found that the results of MA05 models imply younger ages 
by a factors up to 6 and lower stellar masses, by 60% on aver- 
age, with respect to BC03 (and in general to the others models of 
population synthesis). The overestimate of BC03 galaxy stellar 
masses, produced by the lack of the TP-AGB contribution to the 
integrated luminosity, become considerable for evolved stellar 
populations with ages between 0.5 and 2 Gyr. 

In Fig. 6 we compare the stellar masses estimated for the 
z > 3.5 galaxies with the two different template libraries (BC03 
vs MA05). It appears that the BC03 models tend to overestimate 
galaxy masses of our sample with respect to the MA05 ones, by 
0.15 dex on average (see the smaller panel on the top of Fig. 6), 
in good agreement with the previous studies (e.g. Maraston et al. 
2006; Wuyts et al. 2007; Cimatti et al. 2008; Werner et al. 2004). 
In the rest of this paper we will adopt the MA05 stellar tem- 
plates to compute the galaxy stellar masses and the comoving 
stellar mass density. However, in some cases we also use the 
BC03 models results for comparisons with the literature. 
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Fig. 6. Comparison between estimated masses by means of different 
template libraries for our final sample of 53 objects: log M(BC03) vs log 
M(MA05). It is visible that BC03 models tend to overestimate galaxy 
masses with respect to the MA05 ones. The inserted panel shows the 
distribution of the log M(BC03)-log M(MA05) relation as a function of 
the galaxy mass. 



5.2. Comparison with Lyman break selection techniques 

One of the most popular method used to identify high-z galaxies 
is the 'dropout' technique, a colour selection based on the red- 
shifted 912 A Lyman-break caused by neutral hydrogen in the 
galaxy SEDs. This technique was introduced by Guhathakurta 
et al. (1990) and Steidel & Hamilton (1992) to select Lyman 
Break Galaxies (LBG) at z - 3 (u-band dropout). It has been 
largely used also to identify galaxies in the redshift range z ~ 
4 - 6, as b-, v- or i-dropout (Steidel et al. 1999; Dickinson 1998; 
Steidel et al. 2003; Giavalisco et al. 2004b). However this colour 
criterion is strongly biased against red massive galaxies at high 
redshift, that are faint or totally undetected in optical and near-IR 
bands. 

To obtain a quantitative estimate of the fraction of red ob- 
jects missed by combining optical magnitude limited selection 
and dropout techniques, we matched our final sample of IRAC- 
selected high redshift galaxies (53 objects) with a sample of B- 
and V-dropout LBGs in the GOODS-N field from Giavalisco 
et al. (2004a). We found that only 7 of our massive high-z 
candidates are included in the B-dropout sample, and none in 
the V-dropout sample. This means that the remaining IR bright 
(rri4.5 < 23) z > 3.5 candidates have been missed by those se- 
lection criteria, due to their faint emission at optical and near-IR 
wavelengths (most of them having zsso > 27). 

While the vast majority of our candidates are missed by the 
dropout selection technique, there could also be a significant 
population of massive galaxies with m\ .5 < 23 that are identi- 
fied through the dropout technique and that have been missed by 
our selection criteria. We tested for the presence of such galaxies 
using again the sample of Giavalisco et al. (2004a). 

Giavalisco et al. (2004a) listed 684 LBGs candidates (B- 
dropout + V-dropout). Only 7 of those are also part of our sam- 
ple of z > 3.5 candidates. Only 32 of the Giavalisco et al. galax- 
ies have m 4 . 5 < 23 (15 B-dropout and 17 V-dropout). Of these, 
20 were excluded from our sample because of their photomet- 
ric redshifts z < 3.5 (8 B-dropout and 12 V-dropout). One more 
object (V-dropout) was excluded because of the presence of X- 
ray emission and 4 objects ( V-dropouts) were excluded because 
of the lack of an IRAC-8.0 yum peak. One of these 4 V-dropout 
galaxies has a confirmed spectroscopic redshift of z = 5.2 (see 
also Sect. 5). It is reasonable to suppose that all the 4 V-dropouts 
galaxies that we excluded are actually genuine z > 3.5 galaxies 
with young ages, lacking a distinct peak in the IRAC-8.0 //m 
band. Even taking them into account would increase our sample 
of z > 3.5 galaxy candidates by only 7% and would reduce to 
81% the fraction of z > 3.5 massive galaxies that are lost by the 
UV dropout selection techniques. 

It is also interesting to make a comparison with the Lyman- 
break galaxies samples selected by Eyles et al. (2007) and Stark 
et al. (2007) in the GOODS-South field. The former work re- 
ports 6 galaxies at z ~ 6 with masses in the range 1 - 2.4 x 10 10 
M Q , and 4.5 fim magnitude 1114 5 >23. Since the authors used 
Salpeter IMF and BC03 models to derive galaxy SEDs, we have 
to consider that galaxy masses for these 6 objects are overesti- 
mated by a factor of ~ 2.4 (2.4=1.4 x 1.7, being 1.4 and 1.7 
the conversion factors from BC03 to MA05 and from Salpeter 
to Chabrier IMF, respectively Pozzetti et al. 2007; Cimatti et al. 
2008, see also Sect. 5.1 and Sect. 6.2). When rescaled, masses 
span the range 0.4 - 1 x 10 10 M Q , with only one galaxy with 
M~ 10 10 M Q . Lyman Break Galaxies similar to the ones found 
by Eyles et al. (2007) would not be included in our sample, since 
these are fainter than our magnitude cut 1114.5 = 23. This issue is 
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in agreement with the lower masses found for these LBGs with 
respect to our high-z sample. 

The latter work also reports 6 galaxies, spectroscopically 
confirmed at z ~ 4.4 - 5.6. We used the same conversion factor 
to rescale galaxy masses from BC03 and Salpeter IMF and we 
found that only three galaxies of the Stark et al. (2007) sample 
have masses > 10 10 M . Moreover only two out of these three 
sources have 1114.5 < 23. 

Hence only 4 galaxies out of the 12 LBGs at 4.4 < z < 6 
from the two works, fall within our mass bracket (M~ 10 10 - 
10 11 M Q ), and only two would meet our magnitude selection 
(m4.5 < 23). The two IRAC brightest and most massive sources 
from Stark et al. (2007) have ages ~ 150 Myr, if rescaled by a 
factor of 6 to account for the overestimate produced by BC03 
models (Maraston et al. 2006), and E B -v ~ 0.0. We can sup- 
pose that they are similar to the youngest unreddened popula- 
tion excluded from our final sample by the 8 yum peak criterion 
(see Sect. 5 and Fig. 5, top panel), such as the 4 V-band dropout 
galaxies of the sample of Giavalisco et al. (2004b). The other 
two galaxies with masses ~ 10 10 M G and fainter than our IRAC 
magnitude selection would be even younger (ages < 100 Mry). 

In conclusion, it seems that we are not missing a significant 
population of sources with our selection criteria. A few blue star 
forming galaxies with young ages are indeed excluded but these 
objects have generally low stellar masses that would not meet 
our stellar mass completeness limits (M* > 5 x 10 10 M G ; see 
Sect. 6.1 and 6.2). 

5.3. 24 nm emission 

We used observations at 24 fan with Spitzer +MIPS to give some 
constraints on the "activity" of our high redshift galaxy candi- 
dates. Detection of 24 fim flux for z > 3.5 galaxy candidates 
could be explained in terms of radiation coming from star for- 
mation activity lately re-processed and re-emitted by the dust 
(emission lines from PAH molecules), or with the presence of 
an highly obscured AGN in relatively quiescent galaxies. We di- 
vided our sample of high redshift candidates in two sub-samples: 
MlPS-detected (MlPS-d) and MIPS -undetected (MIPS-u) sam- 
ples and discuss them separately in the following. 

5.3.1. The MIPS-u sample 

21 objects in our sample are undetected in MIPS down to the 
limit of 20 /iJy. The best fit Spectral Energy Distributions (SEDs) 
of the MIPS-u sub-sample is presented on the left panels of 
Fig. 13. The right panels of the same figure show the^J distri- 
bution as a function of redshift for each of the sources. We also 
report the most probable photometric redshift solution (z p hot)- 
The best fitting parameters and the estimated galaxy masses are 
presented in Table 3. Most of the MIPS-u redshifts are spread 
around a median value of z ~ 4, in a redshift range z=3.5-5. We 
find that the best fitting SEDs givexl < 2 for 20 out of 21 ob- 
jects and no degeneracy on the redshift solution for the majority 
of the sources. On the contrary, we find also lower redshift solu- 
tions within the 90% confidence intervals of xl distribution for 
5 objects (# 4008, # 9561, # 12327 and #14130), caused mainly 
by the redshift-reddening degeneracy (see Tab. 3 and Fig. 13). 

To further constraints on the nature of the stellar populations 
of MIPS-u galaxies we compared their z - OT4.5 colours with the 
expected ones for MA05 SSP models with different ages (0.1, 
0.3 and 0.5 Gyr). The results are shown in Fig. 7. Although the 
MIPS-u objects span a large range of z - niA,.i colours, most of 



them show colours consistent with older starlight (0.3-0.5 Gyr). 
On the contrary, a smaller part of the sub-sample (7 objects) 
shows bluer colours, more compatible with younger SSPs (0.1- 
0.3 Gyr), probably implying ongoing star formation activity. 

The SED fitting analysis and the lack of mid-IR 24 fim emis- 
sion could thus suggest that many of the MIPS-u galaxies might 
be 'quiescent', non dusty candidates at z > 3.5. However, we 
should not forget that the MIPS observations in the field are not 
deep enough to unequivocally constrain the lack of star forma- 
tion in these galaxies. The 24 fim emission is generally used 
as SFR indicator in high-redshift galaxies (z ~ 2 to 3; Daddi 
et al. 2005; Papovich et al. 2007; Yan et al. 2007; Rigby et al. 
2008). Following Chary & Elbaz (2001) we can estimate that a 
source at z ~ 4 with a 24 fim flux of 20 fiSy will have a total 
IR luminosity of ~ 1.5 x 10 13 L , corresponding to a SFR of 
~ 1500 Moyr 1 for a Chabrier IMF (Kennicutt 1998). Therefore 
the 24 fim non-detection imposes only a very loose upper limit to 
the SFR. However, the 24 fim non-detection is a necessary con- 
dition for quiescent systems (given the MIPS flux limit). In the 
following analysis we will use the MIPS-u galaxy sample to de- 
rive an upper limit to the space density and stellar mass density 
of 'quiescent' candidates at z ~ 4 (Sect. 6.3). 
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Fig. 7. Comparison between the observed z - 1114 5 colours of the MIPS- 
u galaxies (black open diamonds) and the ones predicted at the same 
redshift for SSPs models (MA05) with different ages (black curves in 
the figure, with age of 0.1, 0.3 and 0.5 Gyr, respectively, from the lower 
part upward). Most MIPS-u objects have z-m 4 S colours consistent with 
SSPs with ages in the range of 0.3-0.5 Gyr. On the contrary, a smaller 
part of the sub-sample (7 objects) shows bluer colors, more compatible 
with younger SSPs (0.1-0.3 Gyr), probably exhibiting star formation 
activity. 



5.3.2. The MlPS-d sample 

Out of the 53 galaxies of our high-z sample, more than half 
(32) are detected at 24 fim. The best fitting SED results for the 
MlPS-d objects and their relative^ distributions within the 90% 
confidence intervals are shown in Fig. 14. In Table 4 we report 
the best fitting parameters for the MlPS-d sample. The quality 
of the fit for the MlPS-d SEDs is very good, with xl values 
below 1.5 for almost all objects (except for #5511, which has 
xl = 2.46). However most galaxies in the sub-sample present de- 
generate redshift solutions, caused by the different models used 
and the amount of reddening present (# 702, # 2560, # 2796, 
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#ID 


Zphot 


A v 


P 2 

xi 


V 


<-inf— 90% 


Zsup— 90% 


Zl-phot 


P 2, 


log M(MA05) 


#522 


3.928 


0.410 


94.28 


0.80 


3.060 


5.230 


10.060 


43.52 


9.960 


#1098 


4.369 


1.465 


14.55 


0.20 


4.250 


4.460 


4.922 


0.04 


10.143 


#251 A 


3.620 


0.307 


97.98 


0.40 


2.570 


4.110 


6.098 


0.47 


10.581 


#3081 


3.543 


0.466 


91.26 


0.70 


3.340 


3.760 


0.456 


0.00 


9.816 


#3302 


4.292 


0.306 


98.01 


0.30 


3.620 


4.600 


10.060 


0.00 


10.565 


ii a r\r\a 

#4008 


9.024 


0.201 


99.63 


0.40 


2.150 


10.060 


3.753 


99.57 


1 r\ r\f (\'A 

10.069 


#5073 


4.481 


0.738 


68.90 


0.20 


4.250 


4.670 


3.417 


20.62 


11.073 


#5460 


4.236 


0.854 


57.65 


0.00 


3.480 


4.390 


3.683 


50.87 


10.402 


#6876 


3.620 


0.683 


74.09 


0.20 


2.220 


4.040 


3.172 


72.67 


10.485 


#7042 


3.774 


0.646 


77.54 


0.50 


3.340 


3.900 


3.529 


75.10 


9.723 


#7286 


3.564 


5.582 


0.00 


0.10 


2.850 


3.620 


3.193 


0.00 


10.610 


#7309 


3.613 


1.267 


24.30 


0.00 


3.130 


3.830 


3.417 


22.90 


9.723 


#7785 


3.970 


1.964 


3.28 


0.30 


3.900 


4.040 


0.540 


0.00 


10.391 


#8403 


4.264 


2.018 


2.76 


0.00 


4.110 


4.390 


3.347 


0.57 


10.710 


#9537 


4.103 


0.150 


99.89 


0.40 


3.200 


4.390 


6.651 


0.00 


10.801 


#9561 


4.026 


0.838 


59.14 


0.00 


2.080 


4.320 


3.396 


51.28 


10.412 


#11682 


4.320 


0.367 


96.11 


0.00 


4.040 


4.530 


2.941 


64.89 


10.186 


#12327 


4.971 


0.202 


99.41 


0.70 


2.920 


7.750 


3.256 


96.17 


11.185 


#14130 


4.432 


0.315 


97.06 


0.10 


1.590 


5.090 


3.067 


96.30 


10.235 


#14260 


4.152 


1.243 


25.71 


0.50 


4.040 


4.250 


0.652 


0.01 


9.935 


#15268 


4.775 


2.534 


0.47 


0.80 


4.390 


5.090 


7.260 


0.07 


10.373 



Table 3. Best fitting parameters for the MIPS-u sample. In the order: identification number, z (best fit solution), probability relative to z, 
reddening, upper and lower bound for the 90% confidence interval, z 2 -p(io<(second best solution), probability relative to Zi-phot- We should point 
out that, for galaxies with z p hot ~ 8 — 10, the best fit redshifts result uncertain, and the equally probable solution range (z,„/ - z sup ), within the 90% 
of confidence, is very large. Also for some of these objects -marked by the footnote"- the hyperzmass procedure was not able to find any physical 
solution by fitting the data points with MA05 models the extreme redshifts. 

a This value for the galaxy mass was found for the secondary solution of photometric redshift, z = 3.753, because the hyperzmass software did not find a plausible solution at z = 9.024. 



# 3393, # 3850, # 5 157, # 5367, # 6099, # 9301, # 9306, #10176, 

# 11304, # 13129, # 14105). 

It is difficult to understand the nature of these objects, given 
that often high and lower photometric redshift solutions result 
equally probable. In these cases the 24 /itn emission could be 
interpreted in different ways, depending on the redshift of the 
objects and on the reddening parameter. Assuming that z > 3.5 is 
correct, and taking into account the lack of hard-X ray emission, 
the MIPS flux could ensue from highly obscured AGNs hosted in 
relatively quiescent galaxies. However, still consistently with the 
high redshift solution, the 24/im emission could be also due to 
star formation at extreme high rates. In the next section we will 
examine this possibility using sub-mm/mm photometry from the 
SCUBA 'supermap'(Borys et al. 2003; Pope et al. 2005), and the 
MAMBO (Greve et al. 2008) and AzTEC (Perera et al. 2008) 
imaging surveys of GOODS-N. 

On the other hand if the lower redshift solution (z ~ 2 - 3) is 
considered, the 24 ^m emission could be explained in terms of 
PAH emission from dust heated by star formation activity, and 
those objects would be accounted as dusty star-burst galaxies at 
z ~ 2 - 3. 

5.4. Sub-mm detections : a large population of massive 
starburst galaxies atz ~ 4? 

Sub-mm/mm selected galaxies (hereafter SMGs) are the bright- 
est star forming galaxies known, being much more more lumi- 
nous than the local ultra-luminous IR galaxies -ULIRGs- with 
Ljr » 10 12 L Q . These galaxies are massive young objects seen 
during their formation epochs, with very high SFR (Lilly et al. 
1999; Scott et al. 2002), about one order of magnitude larger than 
that of typical systems with similar masses (Daddi et al. 2007). 
They are also fairly rare objects, probably due to the short du- 



ration of their bright phase (« 100 Myr; Greve et al. 2005). 
These objects could represent the common early phase in the 
formation of massive elliptical galaxies, and hence they might 
be crucial link to understand the massive galaxy formation pro- 
cess. 

The numerous sub-mm datasets available in the GOODS-N 
field can be used to improve the constraints on the nature of our 
z > 3.5 candidates. 

We made direct comparison with three sub-mm maps of the 
field: 

- the Sub-millimeter Common User Bolometer Array 
(SCUBA, Holland et al. 1999) 'supermap' of 35 SMGs 
selected with S/N> 4 at 850 fim (Borys et al. 2003; Pope 
et al. 2005, 2006); 

- the Max Planck Millimeter Bolometer Array (MAMBO, 
IRAM 30-m telescope, Bertoldi et al. 2000) 1200 //m map 
from Greve et al. (2008), composed by 27 sources detected 
with 3.5 < S/N < 4 and 20 with S/N > 4; 

- the AzTEC (Wilson et al. 2008, ,15 m James Clerk Maxwell 
Telescope -JCMT) 1 100 /im map of 28 sources detected with 
S/N < 3.75 (Perera et al. 2008). 

We used a simple approach, looking for overlap of our sam- 
ple with the position of the submm/mm selected galaxies in the 
above works. We have used radial separation limits of 7", 6" and 
9", to search for matches with galaxies in the SCUBA, MAMBO 
and AzTEC, respectively. 

We found eight of our high-z candidates in the error box of 
the SCUBA 'supermap' sources (#5157, #6099, #6463, #8520, 
#10176, #11682, #13857, #15541), but only seven (except 
#1 1682) correspond to the IRAC counterparts identified by Pope 
et al. (2006). For two of these galaxies, GN20 and GN20.2, 
(objects #15541, #13857 in our sample) spectroscopic redshifts 
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y z 

A y 


P 2 


A v 


mf- 


SUp- c_ 


<-2-phot 


P !, 


log M(MA05) 


# 588 


4.285 


0.785 


64.36 


0.60 


3.690 


4.740 


9.388 


0.05 


9.973 


# 702 


8.814 


0.067 


99.99 


0.40 


1.450 


10.060 


3.340 


99.99 


10.563 


#1720 


3.522 


0.526 


87.33 


0.60 


2.850 


4.250 


3.921 


84.30 


10.399 


#2172 


3.711 


0.536 


86.56 


0.70 


3.550 


4.180 


4.117 


70.42 


1 1 .097 


#2560 


3.830 


0.044 


100.00 


0.40 


2.780 


8.520 


7.274 


100.00 


10.762 


#2796 


10.060 


0.279 


98.59 


0.80 


3.360 


10.060 


3.536 


95.36 


1 r\ o'irx',} 

10.870 


#2894 


3.893 


0.628 


79.14 


0.70 


2.990 


4.460 


10.060 


0.37 


10.441 


#3393 


5.503 


0.477 


90.59 


0.80 


3.830 


5.860 


5.167 


84.73 


10.647 


#3850 


4.166 


0.063 


100.00 


0.70 


3.410 


9.290 


8.009 


99.90 


10.940 


#5157 


4.740 


0.344 


96.90 


0.00 


4.110 


7.400 


6.504 


92.59 


11.433 


#5367 


4.355 


0.494 


89.51 


0.80 


3.900 


4.880 


10.060 


0.03 


10.806 


#5511 


4.152 


2.521 


0.50 


0.50 


4.040 


4.250 


1.310 


0.00 


10.551 


#6099 


5.251 


0.367 


96.10 


0.10 


3.200 


8.450 


4.999 


96.02 


1 1.178 


#6245 


5.062 


0.108 


99.98 


0.00 


4.460 


5.370 


2.934 


81.31 


10.662 


#6463 


9.437 


0.193 


99.68 


0.10 


3.690 


10.060 


4.215 


99.64 


9.973" 


#8520 


3.683 


0.573 


83.71 


0.00 


3.340 


3.830 


3.137 


44.99 


10.730 


#9007 


3.669 


1.117 


34.42 


0.60 


2.780 


4.110 


2.969 


34.41 


10.824 


#9301 


8.226 


0.059 


100.00 


0.20 


2.640 


10.060 


8.037 


100.00 


11.035 


#9306 


8.660 


0.127 


99.95 


0.80 


2.360 


10.060 


8.219 


99.92 


11.011 


#9782 


3.928 


1.472 


14.28 


0.20 


3.480 


4.320 


8.723 


0.02 


11.116 


#10176 


5.818 


0.002 


100.00 


0.20 


3.550 


9.570 


5.517 


100.00 


1 i ion 

1 1.280 


#1 1304 


4.145 


0.378 


95.67 


0.80 


3.550 


8.100 


5.461 


89.78 


10.759 


#12699 


3.543 


1.760 


6.21 


0.30 


3.200 


3.690 


3.354 


5.62 


10.558 


#13129 


7.960 


0.030 


100.00 


0.20 


3.410 


10.060 


8.296 


100.00 


1 1.056 


1 1 1 ion 

#13857 


3.746 


0.863 


56.74 


0.10 


3.270 


4.040 


3.403 


43.23 


10.707 


#14505 


3.564 


0.014 


100.00 


0.80 


2.920 


10.060 


7.757 


100.00 


10.025 


#14668 


4.600 


0.209 


99.56 


0.10 


1.870 


5.230 


4.789 


99.53 


10.177 


#14722 


4.145 


1.816 


5.23 


0.40 


4.110 


4.180 


0.617 


0.36 


10.180 


#14793 


4.474 


1.225 


27.41 


0.00 


4.040 


5.160 


9.185 


0.01 


11.411 


#15541 


3.949 


1.216 


27.45 


0.40 


3.690 


4.110 


0.680 


0.00 


11.146 


#15761 


4.292 


0.177 


99.78 


0.10 


3.900 


4.530 


6.567 


0.00 


11.126 


#15771 


3.564 


0.360 


96.35 


0.20 


3.060 


3.830 


6.854 


0.00 


11.283 



Table 4. Best fitting parameters for the MlPS-d sample. In the order: identification number, z (best fit solution), chi 2y , probability relative to z, 
reddening, upper and lower bound for the 90% confidence interval, zi- p hot( second best solution), probability relative to zi-phot- For the objects with 
extremely high redshift best fit solutions, the same caveat we gave in Tab. 3 should be taken into account 



a These values for the galaxy masses were found for the secondary solutions of photometric redshift, because the hyperzmass software did not find a acceptable solution at z ~ 9 - 10. 



around z ~ 4 were recently determined by Daddi et al. (2008, 
hereafter D08). Those objects are among the most distant SMGs 
spectroscopically confirmed known at the present. Our photo- 
metric redshifts for these galaxies are consistent within the error 
(~ 0.2) with the spectroscopic ones (being 3.75 and 3.95 for 
GN20 and GN20.2 respectively). It is interesting to remark that 
D08 discuss the evidence of a galaxy proto-cluster at z ~ 4 in the 
GN20 and GN20.2 area, due to the finding of an over-density of 
B-dropout galaxies. Four objects in total from our z > 3.5 sample 
(#13857, #14722, #15268, #15541) are also within a 25" radius 
centered on the coordinates of GN20. This represents an over- 
density of a factor of 18 with respect to the expected number 
density of IRAC-selected galaxies in the field. Massive galaxies 
seem to be tracing this proto-cluster structure at z = 4.05. 



We found five sources of our high-z sample matching the 
MAMBO galaxies within the correlation distance. Out of them 
only #15541 (GN20) is also SCUBA-detected. In addition, 8 
objects from the AzTEC survey correspond to galaxies in our 
Z > 3.5 sample, within the adopted search radius. Three of them 
are also SCUBA detected, and three are in common with the 
MAMBO survey from G08. 



Candidate z > 3.5 galaxies in our sample that are consis- 
tent with being counterpart of SCUBA, MAMBO and/or AzTEC 
sources are listed in Tab. 5. 

It should be noted that the bulk of the sub-mm counter- 
parts, except the objects #522 (GN1200.5/AzGN06) and #7785 
(GN1200.46), are 24 yum-detected (MlPS-d sample). The latter is 
also a less secure MAMBO detection, being among the sources 
with a S/N below the 4.5 cr limit, and is not part of the 'robust 
deboosted catalogue' . 

We can check the reliability of the high-z solution for 
the sub-mm detections by comparing their S850/S24 vs S24/S8.0 
colours with the ones found by D08 for the two confirmed SMGs 
galaxies at redshift ~ 4 (GN20 and GN20.2). We show this com- 
parison in Fig. 8, where GN20 and GN20.2a are represented 
as black bold diamonds. As suggested by D08, starburst galax- 
ies at z ~ 4 would have relatively blue Spitzer MIPS-IRAC 
colours, similarly to GN20 and GN20.2a (log(S 2 4/S 8 .o) < 0.7, 
log(S24/S4.5) < 1.11). Moreover, because of the negative sub- 
mm k-correction the S850/S24 ratio should increase with redshift. 
Hence we expect that starburst at z ~ 4 would occupy the upper- 
left part of the diagram in Fig. 8, with log(S85o/S24)> 2 (such as 
GN20 and GN20.2a). 
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4.14 


3.55 
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2.2 


0.8 


#6099 
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5.25 


2.95 


10.0 


2.0 


1.9 


1.8 




2.5 


#10176 
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.../AzGN31 


5.818 


3.55 


10.0 


2.4 


13.42 


8.9 


2.13 


3.9/5.3 


#5157 
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4.74 


4.11 


7.40 


2.3 


2.45 
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2.4 
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3.34 
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3.20 


8.6 


3.83 
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2.50 


3.2 




3.2 
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GN1200.1/AzGN01 
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3.69 


4.11 


2.95 


4.055 


20.3 


9.3/10.69 


1.3/1.9/1.4 
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3.27 


4.04 
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Table 5. SCUBA, MAMBO and AzTEC counterparts found in our z > 3.5 of MlPS-d detected sources. In the columns we report: the ID number 
of the objects (this paper); the names from the sub-mm surveys; photometric redshift estimated in this paper and their 90% confidence intervals; 
photometric redshifts by Pope et al. (2006, P06); radio-IR redshifts by Daddi et al. (2008, D08) - with the two spectroscopic values in boldface- 
; the SCUBA, MAMBO and AzTEC fluxes. The last column is the distances between the proposed IRAC-4.5 yum counterparts and the sub-mm 
positions. 



a Not deboosted flux 



Out of the SCUBA sources (black diamonds), two 
(#10176/GN09 and #8520/GN12) have mid-IR colours compa- 
rable with GN20 and GN20.2a, and log(S 8 5o/S 2 4)> 2). This find- 
ing can be seen as a further confirmation of the z > 3.5 redshift 
estimated for these objects. It is also in agreement, within the 
99% confidence interval, with the radio-IR redshifts estimated 
by D08 (see their Tab. 3). In the same figure also the 'SCUBA- 
blank' MAMBO (red triangles) and AzTEC (green crosses) de- 
tections are shown. To represent them in this diagram we scaled 
the S1200 and Si 100 fluxes to Ssso, multiplying by a factor of 2, i.e. 
the typical S 85 o/Si2oo ratio found for GN20 and GN20.2 (D08). 
This is also consistent with the Ssso/Snoo flux ratio derived by 
Perera et al. (2008, 2.08 + 0.18). For the MIPS-u #522 and #7785 
we show the limits in the two-colours diagram, computed with 
S 2 4=20 //Jy (the 5 cr 24 yum-MIPS flux lower limit). Also the 
sub-mm-IR colours of these sources are consistent, within the 
error bars, with the expected colours for high-z SMGs. For com- 
parison we show also the positions occupied in Fig 8 by the 
other sub-mm galaxies from the (Pope et al. 2006) sample with 
secure IRAC counterparts, that do not overlap with our sam- 
ple. Almost all these sources (violet small squares) are spec- 
troscopically confirmed at z ~ 1 - 2.5 (see Pope et al. 2006). 
The evidence that most of them occupy the lower-right part of 
Fig 8 strongly support the above statements. Other two galax- 
ies from the SCUBA 'supermap' have colours similar to GN20 
and GN20.2, i.e. GN10, and GN22. The first one is claimed by 
several authors (D08, Dannerbauer et al. 2008; Wang et al. 2007, 
2009) to be a z > 3.7 starburst, although (Pope et al. 2006) report 
a z = 2.2 photometric redshift. The second one is a spectroscopi- 
cally confirmed galaxy at z = 2.509 (Chapman et al. 2005; Pope 
et al. 2006). It represents the only case in the SCUBA sample 
of a lower-redshift sub-mm/mm emitter with both S850/S24 ratio 
and IRAC colours similar to spectroscopic confirmed z > 3.5 
SMGs. Nevertheless from a statistical point of view the diagram 
shown in Fig 8 seems to be a valid diagnostic to prove high-z 
starburst for sub-mm/mm sources. 

Hence we can conclude that, for ~ 78% of the galaxies 
in our sample that are likely sub-mm/mm detected, these re- 
sults support the hypothesis of extreme SFR activity at high-z. 



Instead the sub-mm (SCUBA) sources lying in the right part of 
Fig. 8(objects #5157, #6099 and #6463 — 22% of the sub-mm 
detections) are more probably at lower redshifts, as proposed by 
Pope et al. (2006). 
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Fig. 8. S850/S24 vs S24/S8.0 colours for the SCUBA(black diamonds), 
MAMBO(red triangles) and AzTEC(green crosses) detections of our 
high-z sample. The comparison with the GN20 and GN20.2 (Daddi 
et al. 2008) seems to support the high-z hypothesis for most of these 
sources (left part of the diagram; see Sect.5.4). The solid black line rep- 
resent the colour limit suggested by (Daddi et al. 2008) for sub-mm 
galaxies at z > 4 (log(S94/S8.o) < 0.7). The small violet squares are the 
SCUBA 'supermap' sources that do not overlap with our high-z sample, 
most of them spectroscopically confirmed at z ~ 1-2.5. 



5.5. Mid-IR vs IR colours 

The properties and the overall spectral energy distributions make 
the MlPS-d galaxies comparable to one of the most debated 
objects in the literature, HUDF-JD2. This object, identified as 
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Fig. 9. Infrared colours mg.o - "J24 vs /W3.6 - mg.o of our high-z sam- 
ple of galaxies (black diamonds) compared with the objects found by 
Rodighiero et al. (2007) (red open triangles). The bold red triangle rep- 
resents HUDF-JD2 (Mobasher et al. 2005). Bold black open diamonds 
are the SCUBA detected galaxies of our sample. The two objects with 
confirmed z ~ 4, GN20 and GN20.2, are highlighted (filled black dia- 
monds). The MlPS-undetected objects are shown as upper limits in the 
mg.o _ m 2A colour: black and red arrows for our sample and Rodighiero 
et al. (2007) sample, respectively. Two evolutionary colour tracks of an 
ULIRG (violet dot-dashed curve; Chary & Elbaz 2001), and of a galaxy 
hosting an obscured AGN (dashed green curve; Seyfert 2 galaxy NGC 
1068), with increasing redshift values, are also shown. The black hori- 
zontal line represents the upper limit suggested by D08 and also shown 
in Fig. 8 for colours of starburst galaxies at z > 3.5. 



a z ~ 3.4 galaxy by Yan et al. (2004) and Chen & Marzke 
(2004), was recently suggested by Mobasher et al. (2005) as 
a massive post star-burst galaxy at z - 6.5 hosting an ob- 
scured AGN. However subsequent studies (Dunlop et al. 2007; 
Rodighiero et al. 2007; Chary et al. 2007) claimed that HUDF- 
JD2 is more likely a lower redshift (z ~ 1.5-2.2) star forming 
galaxy with very strong dust reddening. Our galaxies are also 
similar to candidate z > 4 galaxies found by Rodighiero et al. 
(2007) and Wiklind et al. (2008). In Fig. 9, the infrared colours 
(ms.o - vs /M3.6 - m%.o) of our high-z sample (open black 
diamonds) are compared with the objects found by Rodighiero 
et al. (2007, open red triangles). The bold red triangle represents 
HUDF-JD2 (Mobasher et al. 2005). Bold black open diamonds 
are the sub-mm detected galaxies of our high-z sample. The two 
objects with confirmed z ~ 4, GN20 and GN20.2, are repre- 
sented by filled black diamonds. The MIPS -undetected objects 
are shown as upper limits in the m%x> - '«24 colour (black and 
red arrows for our sample and the sample of Rodighiero et al. 
2007, respectively). The evolutionary colour tracks of an ULIRG 
(Chary & Elbaz 2001) and of a galaxy hosting an obscured AGN 
(the Seyfert 2 galaxy NGC 1068) with increasing redshift values 
are also shown. It is clear, though, that such locally calibrated 
templates are not very useful to interpret the colors of distant 
galaxies. This can be argued just considering that the observed 
ms.o -m.2A colours of GN20 and GN20.2 are ~ 1 mag redder with 
respect to the ULIRG track at z ~ 4. As discussed in Daddi et 
al. (2008), this is most likely due to a lower stellar mass to light 
ratio in the distant galaxies. 

The black horizontal line represents the upper limit sug- 
gested by D08 and also shown in Fig. 8 for colours of star- 
burst galaxies at z > 3.5 (log(S24/S8.o) < 0.7, in AB magnitude: 
W8.o - '«24 < 1-7; see also Sect. 5.4). If the high-z solution is 



correct, objects bluer than this limit should be strongly domi- 
nated by star formation, such as GN20 and GN20.2. The issue 
that most of the the sub-mm detected galaxies in our sample are 
down to the black line in the figure also suggest that they have 
properties similar to the two starbursts spectroscopically con- 
firmed at z ~ 4. The MlPS-d sources which have similar colours 
but are lacking in a sub-mm counterpart could be also high-z 
starbursts. Alternatively, the lack of sub-mm detection should be 
justified either by the shallow sensitivity of the currently avail- 
able sub-mm instruments (SCUBA, MAMBO, AzTEC), or by 
particularly "warm" SEDs for this objects. 

On the other hand, the z > 3.5 candidates with _ m 2A > 
1.7 in Fig. 9 could represent a different population. Their red- 
der colours, more similar to the bulk of the Rodighiero et al. 
(2007) sample, could be due to heavily obscured AGN, as also 
suggested by Rodighiero et al. (2007) for most of their sources, 
or they might be z < 3.5 contaminants. The position of object 
HUDF-JD2, overlapping with the solid black line in the middle 
part of the diagram, does not allow us to give any further con- 
straint on its nature. 

6. The comoving stellar mass density at 3.5 < z < 5 

Recent studies (Drory et al. 2005; Fontana et al. 2006; Stark et al. 
2007; Verma et al. 2007; Eyles et al. 2007; Yan et al. 2006) have 
extended measurements of the total stellar mass density up to 
z ~ 4 - 5, using deep multicolour data from optical and near-IR 
selected samples. One of the goals of this paper is to determine 
how our IRAC-4.5 selected sample contributes to the comoving 
stellar mass density at z > 3.5. In fact, as already discussed, 
observations of the high redshift Universe are generally biased 
against red galaxies with large mass-to-light ratios, potentially 
missing galaxies that could potentially make a strong contribu- 
tion to the total mass density (see Fontana et al. 2006). Here, we 
emphasize that we can only determine a lower limit to the co- 
moving stellar mass density in the redshift interval 3.5 < z < 5, 
for several reasons. First, the combination of our IRAC 4.5 fim 
magnitude limit and the 1.6 //m peak SED criterion leads to 
missing objects with lower masses at these redshifts, such as 
the fainter, bluer Lyman break galaxies, as we have discussed 
in Sect. 5. More generally, a magnitude-limited sample (selected 
at any wavelength) is not the same as a mass-limited sample, as 
we will discuss further in the next section. 



6.1. Incompleteness effects and mass selection criteria 

In estimating the comoving stellar mass density we must take 
into account that our sample could suffer from incompleteness. 
As several other authors (Fontana et al. 2004; Drory et al. 2004; 
Labbe et al. 2005) have noted, IR-selected samples are still not 
equivalent to mass-selected samples. In fact, at any redshift, the 
galaxies detected above the sample magnitude limit can have a 
fairly large range of possible M*/L ratios, depending on their 
spectral properties, ages, dust extinction, metallicities, etc. The 
effect is that magnitude-limited samples, at higher redshifts, are 
progressively biased against objects with lower masses and high 
M*/L ratios, such as galaxies that are not currently forming 
stars, or which are highly extincted. 

We have attempted to evaluate the consequences and mini- 
mize the effects of incompleteness by determining the threshold 
mass limit as a function of redshift for a variety of galaxy mod- 
els spanning a range of ages and degrees of extinction. For this 
purpose we used the MA05 template library with the Chabrier 
IMF, the same as was used in the preceding sections to derive 
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the galaxy stellar masses. We consider both the simplest case of 
SSP models with different ages (0.1, 0.5 and 1 Gyr), and a set 
of models with a constant star formation rate, an age of 0.5 Gyr, 
and different amounts of extinction {Eb-v = 0.0, 0.3, 0.5, applied 
using the Calzetti attenuation law; Calzetti et al. 2000). The ob- 
served 4.5 fim flux as a function of redshift was derived from 
these models, and from this we translate our 4.5 fim magnitude 
limit (»J4 5 < 23) into a limiting stellar mass at each redshift. 
The results are shown in Fig. 10. For passively evolving galaxies 
(i.e., the SSPs) with ages less than 0.5 Gyr, and for star-forming 
galaxies with Eb-v < 0.3, our sample should be fairly complete 
at M > 5 x 10 10 M Q out to z < 5. This is another advantage of 
our 4.5pm sample selection from the very deep GOODS IRAC 
data, which is complete to significantly lower mass thresholds 
at these very high redshifts compared to Zf-selected samples(see 
Fontana et al. 2006). 

Assuming that the redshift estimate for the objects in our 
high-z sample are correct, we derive their contribution to the co- 
moving number density and to the total stellar mass density, both 
for the total sample (MIPS-u+MIPS-d) and for the MIPS-u sam- 
ple only, in the redshift interval 3.5 < z < 5 as described in the 
following sections. The comoving volume in this redshift range 
over the total GOODS-N solid angle (~ 165 arcmin 2 , coincident 
with the ACS z-band area) is V com ~ 0.725 X 10 6 Mpc 3 . 
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Fig. 10. The redshift dependence of the stellar mass limit for our IRAC 
magnitude-limited sample (m^ < 23), derived from MA05 stellar pop- 
ulation synthesis models. The red curves represent SSP models with 
different ages (0.1, 0.5 and 1 Gyr). The green curves are models with 
constant SFR, an age of 0.5 Gyr, and dust extinction, with three dif- 
ferent values of the colour excess E B -v (0.0, 0.3 and 0.5). The points 
indicate the objects of our z > 3.5 sample, with error bars indicating the 
68% confidence ranges on their stellar masses. 



Given this potential incompleteness, and the fact that we are 
undoubtedly missing a large percentage of bluer, lower-mass 
LBGs in the same redshift range (as discussed in section 5.2), 
we must consider our estimate for the stellar mass density at 
3.5 < z < 5 to be a lower limit. Fig. 10 shows the derived stel- 
lar masses versus redshift for the objects in our IRAC-selected 
sample. We count objects and sum their stellar masses for galax- 
ies above our mass threshold (M~ 5 x 10 10 M ). The resulting 
comoving number density and the stellar mass density for this 
sample are 2.6 x 10~ 5 Mpc~ 3 and (2.9 ± 1.5) X 10 6 M Mpc~ 3 , 
respectively. 

In Fig. 1 1 we show our lower limit to the comoving stellar 
mass density (p„) at 3.5 < z < 5, and compare it to another 
estimate (Fontana et al. 2006, F06) at 3 < z < 4 for galaxies 
in a similar mass range (M* = 3 x 10 10 - 3 x 10 11 M G ). We 
take the F06 value from the compilation published by Wilkins 
et al. (2008), and have corrected it for the M/L differences be- 
tween the Salpeter and Chabrier IMFs (log M(Chabrier)=log 
M(Salpeter)-0.23), as well as by a factor of 1 .4 to account for the 
difference between the BC03 used by F06 and the MA05 models 
used in our study (Sect. 5.1). We also show a curve derived from 
the Millennium Simulation models for the redshift evolution of 
the stellar mass density, computed for galaxy masses above the 
threshold limit of > 5 x 10 10 M Q . Assuming that the photometric 
redshifts for all the selected high-z candidates are correct, our 
results are in good agreement with the Millennium Simulation 
predictions. 




z 



Fig. 11. Comparison between our estimate of the lower limit of the stel- 
lar mass density (p*) at 3.5 < z < 5 (red triangle) and the stellar mass 
density estimated by Fontana et al. (2006, F06) at 3 < z < 4, for galax- 
ies in the mass range M t ~ 3 x 10 10 - 3 x 10" M . The curve shows 
the evolution of stellar mass density with redshift from the Millennium 
Simulation models, computed for galaxies above our threshold limit of 
> 5 x 10 10 Mo. 



6.2. Contribution to the stellar mass density from our total 
sample in the redshift bin z = 3.5 - 5 

Based on the preceding analysis, we assume a limiting stellar 
mass M~ 5 x 10 10 M as the completeness threshold for our 
magnitude-limited sample. For quiescent galaxies with interme- 
diate ages (~ 0.5 Gyr), we should be complete above this mass 
threshold out to z ~ 5. However, for z > 4, we may be incom- 
plete near this mass limit for the oldest possible passive stellar 
populations, or for heavily dust-obscured galaxies. 



6.3. Stellar mass density for the MIPS-u candidates in the 
redshift bin z = 3.5-5 

As described above, our sample should be largely complete for 
'quiescent', non-star-forming galaxies with ages < 0.5 Gyr and 
masses M > 5 x 10 10 M Q . However, at these large redshifts, 
the flux limit at 24 pm, even for very deep data from GOODS, 
corresponds to quite large star formation rates > 1000 M Q /yr. 
Therefore, we cannot affirm that the MIPS -undetected (MIPS- 
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u) objects in are sample are truly 'quiescent', and thus we can 
only derive an upper limit for the number and stellar mass densi- 
ties of quiescent galaxies. For the MIPS-u subsample with stel- 
lar masses > 5 x 10 10 M at z = 3.5 - 5, we find a comoving 
number density of 0.97 x 10~ 5 Mpc~ 3 , and a stellar mass density 
of (1.15 + 0.7) x 10 6 M Mpc~ 3 (MA05 templates and Chabrier 
IMF). 

In Fig. 12 we compare our results for the MIPS-u sample 
with estimates of the stellar mass density for quiescent galax- 
ies at other redshifts from the literature. To facilitate this com- 
parison, here we have scaled all values to the Chabrier IMF 
and MA05 models by means of the M/L conversion relations 
given above (see Sect.6.2). We have also computed the comov- 
ing stellar mass density of the MIPS-u galaxies considering the 
unlikely case that they were all truly "passive", consistent with 
single burst stellar populations (SSPs) with no dust. The masses 
were derived using the ratio between the observed IRAC-4.5 fim 
flux for each galaxy and the expected flux for an MA05 SSP 
model at the same redshift, assuming an age of ~ 0.5 Gyr and 
mass normalized to 1 M . In this case, the number density of 
MIPS-u candidates above the mass threshold would be larger 
(1.52 x 10~ 5 Mpc~ 3 ). The corresponding stellar mass density 
would be (1.46 + 0.7) x 10 6 M Q Mpc~ 3 . This value is also rep- 
resented in Fig. 12 by the red dashed line. From this exercise 
we concluded that even if all the MIPS-u galaxies were consid- 
ered as old SSPs, the number density of 'quiescent' galaxies at 
Z = 3.5-5 would remain very small in comparison with that 
found at z < 2 - 2.5 by previous studies. 

The stellar mass density we found for the MIPS-u galax- 
ies in the redshift range z = 3.5 - 5 is ~ 2 orders of magni- 
tude lower than the local value (black asterisk; Baldry et al. 
2004), and more than 10 times lower than the values found at 
Z ~ 1.7 by Daddi et al. (2005, ; black filled rhombus) for the 
passive-BzK galaxies in the Hubble Ultra Deep Field (HUDF) 
and at z ~ 2.5 by Labbe et al. (2005, black open square) for 
Distant Red Galaxies (DRGs) in the Hubble Deep Field South. 
Moreover our result is also in good agreement with the mass 
density estimated by Wiklind et al. (2008) for z ~ 5.5 Balmer- 
break-selected galaxies without MIPS detections in the GOODS 
South Field at redshift z ~ 5.5 (open rhombus in the figure). 
Considering that our estimate is an upper limit to the comoving 
stellar mass density of massive and quiescent galaxies, we con- 
clude that there is evidence of strong evolution beyond z ~ 2 
(see also Daddi et al. 2005). The decreasing number density of 
quiescent objects at high redshifts could mean that most massive 
galaxies at early epochs were experiencing strong star formation 
activity. 

7. Summary and conclusions 

In recent years, a substantial number of massive galaxies have 
been found at redshifts as large as z ~ 5. However, although the 
range of stellar populations and observed colors among galax- 
ies at all redshifts studied to date has proven to be quite het- 
erogeneous, most surveys at z > 3.5 to date have primarily se- 
lected samples of star-forming blue galaxies, and have been bi- 
ased against red massive objects that are faint or completely un- 
detected in optical and near-IR bands. 

In this work we have studied galaxies selected at IRAC- 
4.5 fim in the GOODS-N field, in a manner that is complemen- 
tary with respect to optical or near-IR selection. This allows us 
to pick out massive, evolved galaxies at high redshift, for which 
the optical rest frame emission (which is most nearly related to 
the stellar mass) is sampled by the redshifted, rest-frame near- 
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Fig. 12. Comparison between the upper limit to stellar mass density (p») 
found in this paper for 'quiescent' galaxies in the redshift bin z = 3.5 - 5 
and results from the literature at other redshifts. All points have been 
scaled to Chabrier IMF and MA05 masses (see Sect. 6.2). Our estimate 
is represented as a red bold triangle. The horizontal red dashed line 
represent the stellar mass density we derive by considering all MIPS- 
u galaxies to be SSPs with age of 0.5 Gyr. The local mass density in 
passive galaxies is shown by the black asterisk (Baldry et al. 2004). The 
black open square is for Distant Red Galaxies (DRGs) at z = 2 - 3 
in the Hubble Deep Field South (Labbe et al. 2005). The black filled 
rhombus is from Daddi et al. (2005) for passive-BzK galaxies in the 
Hubble Ultra Deep Field (HUDF) at z = 1.39-2, and the open rhombus 
is for Balmer-break MIPS undetected galaxies at z = 5-7 in GOODS- 
South from Wiklind et al. (2008). 

infrared portion of the spectrum. On the other hand, this se- 
lection may also be sensitive to 'active' obscured galaxies at 
high redshift, for which IR emission could be reprocessed ra- 
diation from absorbing material surrounding an AGN, or from 
dust heated by star formation, like that seen in dusty starburst 
galaxies such as local ULIRGs. 

Recently, other studies based on near-IR selection have de- 
tected red, massive galaxy candidates at z p hot ^ 4, but have not 
obtained complete samples. In most of these studies, the samples 
have been selected in the /f-band, with the possibility of missing 
TT-undetected objects visible in the longer-wavelength IRAC fil- 
ters (e.g. Mobasher et al. 2005; Wiklind et al. 2008). Rodighiero 
et al. (2007) used IRAC-3.6 fim selection (#23.6 < 23.26) to iden- 
tify a sample of high-z galaxies missed by optical and A'-band 
selection. However they assumed as additional conditions the 
non-detection of their candidates in optical bands, and a detec- 
tion close to the sky threshold limit in the K band (K > 23, AB 
system), which make their sample a priori incomplete. 

Here, we attempt to select a sample of massive high red- 
shift red galaxies that is as complete as possible, by means of 
IRAC selection, that could allow us to recover massive objects 
that may have been missed by previous studies. The most impor- 
tant results of this paper are summarized in the following items: 

• By selecting at IRAC-4.5 fim (m4 5 < 23), we extracted a 
sample of 4 142 objects in the ~ 165 arcmin 2 of the GOODS- 
N field. We performed an SED fitting analysis to estimate 
photometric redshifts for those objects, and found fifty-three 
candidates at z > 3.5, also requiring that peak of the f v stellar 
spectral energy distribution (at 1 .6 fim in the rest-frame) falls 
in the IRAC-8.0 fim band. We excluded unobscured AGN 
from our final sample, identified by their hard and/or soft-X 
ray emission. Almost 8 1 % of our candidates are completely 
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missed by the B- and V-band Lyman break dropout tech- 
niques designed to identify UV-bright, star-forming galaxies 
at similar redshifts. For each object, we evaluate the relia- 
bility and the confidence limits of our photometric redshift 
estimates with a x 1 test. Two objects have spectroscopically 
confirmed redshifts (GN20 and GN20.2; Daddi et al. 2008) 
in good agreement with our estimates. 

• We divide our final sample in two sub-samples of (32) MIPS- 
detected and (21) MIPS -undetected objects (the MlPS-d and 
MIPS-u samples, respectively). 

MlPS-d sample: In the subsample detected at 24/jm, we 
found 18 galaxies with unambiguous high photometric red- 
shift solutions and 14 with degenerate solutions. If the high 
redshift solution is correct, the 24 //m emission could in- 
dicate the presence of a heavily obscured AGN (perhaps 
Compton thick due to the absence of X ray emission), or per- 
haps strong star formation activity in a dusty, hyperluminous 
starburst. We cannot firmly exclude the possibility for these 
galaxies to be dusty starburst galaxies at lower redshifts, for 
which the 24 /im flux could be explained with PAH emission 
atz ~ 2 - 3. 

MIPS-u sample: 21 of our z > 3.5 candidates are unde- 
tected at 24/im. The lack of 24 /im emission is a necessary 
but not sufficient condition to indicate that these are quies- 
cent galaxy. The upper limit to the SFR set by MIPS de- 
tection limit (/24 < 20 /iJy) corresponds to a star formation 
rate < 1500 M Q yr _1 , assuming bolometric corrections de- 
rived from spectral templates for local ultraluminous infrared 
galaxies. For this reason we can use this sub-sample only to 
give an upper limit to the total stellar mass of 'quiescent' 
galaxies at z > 3.5. 

• 14 galaxies from our z > 3.5 sample are detected at submil- 
limeter and millimeter wavelengths with SCUBA (Pope et al. 
2006), MAMBO (Greve et al. 2008) and AzTEC (Perera 
et al. 2008). All but two of are also detected at 24//m. 
Two objects, GN20 and GN20.2, have spectroscopically con- 
firmed redshifts at z ~ 4 Daddi et al. (2008). By comparing 
their colours with other SMGs galaxies in the S850/S24 vs 
S24/S8 diagram, we find that ~ 78% of the sub-mm detected 
objects among our candidates have properties that are fully 
consistent with those of starburst galaxies at z > 3.5. The 
remaining ~ 22% are more likely at lower redshift. 

• We have computed the contribution to the number density 
and to the stellar mass density by the total (MIPS-d+MIPS- 
u) sample and by the subsample that is undetected at 24/im 
(MIPS-u) for galaxies with stellar masses above a limit of 
5 x 10 10 M Q , chosen to minimize incompleteness at z < 5. 
For the redshift range 3.5 < z < 5, we find a comov- 
ing number density 2.6 x 10~ 5 Mpc~ 3 for average masses 
of ~ 10 11 M . The corresponding stellar mass density is 
(2.9 + 1.5) x 10 6 M Q Mpc~ 3 . Our results are in good agree- 
ment with previous determinations in the literature, and with 
the prediction from the Millennium simulations. 

For the MIPS-u sample, we determined a number density of 
0.97 x 10~ 5 Mpc~ 3 and a corresponding stellar mass den- 
sity of (1.15 ± 0.7) x 10 6 M Q Mpc~ 3 . We also considered 
the extreme case that all of the MIPS -undetected galaxies 
were purely quiescent, with ages ~ 0.5 Gyr. This results in 
larger upper limits for their number and stellar mass den- 
sities, 1.52 x 10~ 5 Mpc~ 3 and (1.46 + 0.7) x 10 6 MoMpc" 3 
respectively. Our MIPS-u sample of quiescent candidates ac- 
counts for ~ 5 - 6% of the mass density that was in place at 
z ~ 2 and only for ~ 1% of local stellar mass density (Baldry 
et al. 2004). This provides evidence for strong growth in the 



number of massive galaxies between redshift z = 4 - 5 and 
z = 2. 
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Fig. 13. Spectral Energy Distribution best-fit models for the MIPS-u sample in units of AB magnitude. The red diamonds represent the observed 
magnitudes from the photometric data. Non-detections are represented by upper limits. See also Tab. 3. 
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Fig. 14. Spectral Energy Distributions of MlPS-d galaxies of our final high-z sample (on the left panels), and the correspondent^ distributions as 
a function of redshift (right panels) in 90% confidence intervals. For each candidate the most probable redshifts (first and second redshift solutions 
: Zphot and Zphoti) ar >d the^J values with the respective probability are also labeled. See also Tab. 4. 
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Fig. 15. Multi-wavelength images of the fifty-three galaxies of the final sample (MIPS-d+MIPS-u). 
The cut-outs have a size of 10" x 10". North is up, East is at left. 
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Table 6. Multi-band photometry for the final sample of 44 candidates at z > 3.5. The coordinates reported here are from the Spitzer data-set. 

H Ks mi; 6 "14 5 '«5 8 '«8 '«24 RA Dec 

> 23.0 >23l 23288 2Z895 22~517 2L583 >2(L6 12:36:27.2 62:06:06.1 
>23.() > 23.1 23.379 22.960 22.274 21.836 20.218 12:36:32.7 62:06:21.8 

>23.() > 23.1 22.964 22.857 22.254 20.081 12:36:15.7 62:06:43.4 

23.242 23.127 22.422 22.465 22.065 21.807 > 20.6 12:36:36.8 62:07:14.0 

> 23.0 23.021 22.954 22.611 22.354 22.218 20.161 12:36:51.0 62:08:29.8 
23.187 22.482 21.360 21.118 20.898 20.928 19.156 12:36:13.3 62:09:02.0 

> 23.0 23.222 21.999 21.587 21.227 21.137 18.527 12:35:53.0 62:09:30.2 

> 23.0 23.778 22.958 22.559 22.510 22.183 > 20.6 12:36:04.7 62:09:25.5 

> 23.0 22.404 21.869 21.528 21.167 21.227 18.029 12:36:03.1 62:09:47.5 
>23.0 > 23.1 23.399 22.997 22.608 22.521 20.496 12:37:10.0 62:09:56.4 

>23.0 > 23.1 22.859 22.638 22.480 > 20.6 12:36:44.4 62:10:03.9 

>23.0 > 23.1 22.834 22.638 22.146 22.169 > 20.6 12:35:59.9 62:10:08.4 

23.329 > 23.1 22.423 22.064 21.615 21.376 19.581 12:37:19.5 62:10:21.4 

> 23.0 23.202 22.073 21.579 21.125 21.066 20.037 12:37:12.5 62:10:36.1 

> 23.0 23.330 22.069 21.817 21.405 > 20.6 12:36:53.7 62:11:12.9 

23.058 23.489 22.107 21.831 21.602 21.165 > 20.6 12:36:36.1 62:11:54.7 

> 23.0 23.846 21.479 21.102 20.764 20.856 18.519 12:36:37.5 62:11:57.2 
>23.0 > 23.1 22.729 22.373 21.895 21.745 20.189 12:37:18.9 62:12:17.9 
24.431 24.036 23.019 22.691 22.781 22.567 > 20.6 12:36:31.3 62:12:18.3 
22.865 22.320 21.429 21.444 21.219 20.221 18.084 12:36:37.1 62:12:31.1 
>23.0 > 23.1 22.129 21.678 21.269 21.209 18.305 12:36:08.6 62:12:51.2 
24.283 24.065 22.791 22.688 22.433 22.222 19.809 12:37:20.6 62:13:00.2 

> 23.0 23.287 22.875 22.426 21.948 21.903 19.177 12:37:02.5 62:13:02.6 
>23.0 > 23.1 23.042 22.801 22.468 22.386 > 20.6 12:36:30.4 62:13:33.0 

> 23.0 23.463 23.224 22.992 23.337 23.059 > 20.6 12:37:25.5 62:13:39.8 
24.464 22.810 22.519 22.265 21.941 22.033 > 20.6 12:37:09.5 62:13:50.3 
24.377 23.068 23.031 22.875 23.306 22.830 > 20.6 12:36:41.6 62:13:50.1 
23.265 23.026 22.584 22.616 22.637 22.473 > 20.6 12:36:55.9 62:14:12.5 
22.877 22.841 22.198 21.971 22.106 21.684 > 20.6 12:37:03.8 62:14:51.6 
23.162 22.617 21.867 21.739 21.504 21.488 20.144 12:36:45.1 62:14:48.8 

> 23.0 23.721 22.399 22.033 21.749 21.528 19.433 12:37:34.7 62:15:15.8 
23.802 > 23.1 22.924 22.484 22.332 22.270 20.255 12:38:03.6 62:15:30.5 

> 23.0 22.973 22.360 21.967 21.735 21.555 19.661 12:37:33.8 62:15:32.1 

> 23.0 22.824 22.133 21.873 21.650 21.485 > 20.6 12:36:39.4 62:15:42.8 

> 23.0 23.056 22.686 22.597 22.463 > 20.6 12:37:52.7 62:15:49.3 

> 23.0 23.298 22.203 21.487 21.199 21.046 18.889 12:37:39.5 62:15:59.0 
>23.0 > 23.1 22.428 21.882 21.439 21.237 20.034 12:36:22.0 62:16:16.2 
>23.0 > 23.1 21.918 21.373 20.939 20.617 19.959 12:36:31.9 62:17:14.8 

> 23.0 24.296 22.702 22.682 22.455 22.507 > 20.6 12:37:39.2 62:17:36.8 

> 23.0 24.705 22.645 22.189 22.091 21.367 > 20.6 12:38:05.4 62:18:16.4 
24.640 22.082 22.107 21.959 21.996 21.895 19.942 12:37:24.0 62:18:33.7 
24.368 > 23.1 22.963 22.364 22.029 21.515 20.186 12:37:28.1 62:19:20.5 
24.107 > 23.1 22.605 22.465 21.913 21.578 20.637 12:37:08.8 62:22:02.2 

> 23.0 23.387 23.048 22.900 22.846 22.596 > 20.6 12:37:36.3 62:21:25.3 
24.668 22.678 22.552 22.680 22.814 22.647 > 20.6 12:37:13.0 62:21:11.4 
>23.0 > 23.1 22.924 22.388 22.140 21.573 20.211 12:36:44.0 62:19:38.8 

> 23.0 23.497 22.936 22.749 22.956 22.385 20.610 12:37:12.3 62:23:03.4 
23.167 23.125 22.352 22.714 22.414 19.934 12:37:11.5 62:21:55.8 

> 23.0 23.757 21.952 21.566 20.716 20.492 18.787 12:36:58.9 62:22:15.3 
>23.0 > 23.1 22.790 21.887 21.316 21.071 > 20.6 12:37:09.6 62:22:02.5 

> 23.0 22.753 21.874 21.573 21.010 20.595 19.359 12:37:11.9 62:22:12.4 

> 23.0 22.707 21.688 21.299 20.985 20.931 19.140 12:36:55.3 62:21:07.9 
23.687 22.170 21.418 21.117 20.795 20.743 19.659 12:37:01.5 62:20:25.3 
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